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HYBRID UPWIND SPLITTING (HUS) 
BY A FIELD-BY-FIELD DECOMPOSITION 

Frederic Coquet and Meng-Sing Liou^ 


Abstract 

We introduce and develop in this paper a new approach for upwind biasing : the Hybrid 
Upwind Splitting (HUS) method. This original procedure is based on a suitable hybridiza- 
tion of current prominent Flux Vector Splitting (FVS) and Flux Difference Splitting (FDS) 
methods. The HUS method is designed to naturally combine the respective strengths of 
the above methods while excluding their main deficiencies. Specifically, the HUS strategy 
yields a family of upwind methods that exhibit the robustness of FVS schemes in the cap- 
ture of nonlinear waves and the accuracy of some FDS schemes in the resolution of linear 
waves. 

We give a detailed construction of the HUS methods following a general and systematic 
procedure directly performed at the basic level of the Field by Field (i.e. waves) decom- 
position involved in FDS methods. For such a given decomposition, each field is endowed 
either with FVS or FDS numerical fluxes, depending on the nonlinear nature of the field 
under consideration. Such a design principle is made possible thanks to the introduction 
of a convenient formalism that provides us with a unified framework for upwind methods. 

The HUS methods we propose bring significant improvements over current methods 
in terms of accuracy and robustness. They yield entropy-satisfying approximate solutions 
as they are strongly supported in numerical experiments. Field by field hybrid numerical 
fluxes also achi eve fairly simple and explicit expressions and hence require a computational 
effort between that of the FVS and FDS. 

Several numerical experiments ranging from stiff ID shock-tube to high speed viscous 
flows problems are displayed, intending to illustrate the benefits of the present approach. 
We shall assess in particular the relevance of our HUS schemes to viscous flow calculations. 
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Introduction 


Over the past decades, the research effort in computational fluid dynamics (CFD) 
has resulted in several prominent and widely used schemes. In particular, two distinct 
and complementary approaches for upwind biasing have been proposed and investigated, 
namely the Flux Vector and Flux Difference Splitting methods (FVS and FDS respec- 
tively). Each has demonstrated fairly interesting capabilities, inherited from their design 
principle and assumptions, which are manifested in the numerical calculation in the form of 
either robustness or accuracy. In other words, the design principle directly determines the 
effectiveness of the method; abundant in the literature are papers reporting their successes 
and failures in calculations. 

The first approach (FVS) has received attention mainly for the Euler equations. 
Harten et al. [22] referred it as the Boltzmann approach for its formal links with the 
gas dyn ami c setting where the flow is modeled by a set of pseudo-particles evolving freely 
according to a collisionless Boltzmann equation. As a consequence of the resulting mix- 
ing of particles, FVS methods suffer from a constitutive lack of resolution for linear waves 
[50] (contact discontinuities), thereby making them irrelevant candidates for Navier-Stokes 
calculations since the viscous layer behaves in the limit like a contact discontinuity. On 
the other hand, FVS methods have demonstrated their ability to capture intense and only 
admissible nonlinear waves (such as shock waves). 

The second approach (FDS), generally classified as approximate Riemann solvers, has 
been developed primarily to provide approximate information regarding wave speeds and 
amplitudes in the Riemann problem. Two most outstanding FDS methods are proposed by 
Osher [35] and Roe [45]. These FDS methods specifically provide a direct representation 
of the linear ly degenerate field with which the contact discontinuity is associated. Since 
they both, by construction, reproduce the exact solution for a single stationary contact 
discontinuity, they become accurate and relevant candidates for solving the Navier-Stokes 
(at least the steady) equations. Compared to Godunov’s exact solution of the Riemann 
problem [17], Osher ’s scheme replaces the shock by an overturned compression wave and 
thus makes an explicit computation of numerical fluxes possible, but at the expense of 
reducing the robustness. Jacob [24] reported a lack of robustness in the Osher splitting for 
capturing a strong detached shock. 

More drastic simplification is used in the Roe splitting, rendering it an even cheaper 
procedure in which the transonic expansion is replaced by a shock, but with the conse- 
quence of allowing a physically inadmissible expansion shock. Quirk [41] recently sum- 
marized some failings encounterred by using the Roe splitting in calculating nonlinear 
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waves . 


Therefore, FVS and FDS approaches do not yield the same advantages and drawbacks, 
but they clearly complement each other. 

These observations have prompted some recent works in implementing some switches 
directly aimed at choosing either a FVS or a FDS method according to the local smoothness 
of the solution but at the expense of tunable parameters and of a discontinuous method 
(see [41] for a related idea). 

Up to our knowledge, the first scheme combining both FVS and FDS features into a 
single flux for the sake of robustness and accuracy is the Advection Upstream Splitting 
Method (AUSM) introduced by Liou and Steffen [30]. This original splitting is indeed 
based on the appropriate definition of a cell-face advection Mach number. Liou [32] has 
recently proposed a convenient framework for the development of a family of schemes, 
referred to as AUSM+, which further improves AUSM while preserving its advantageous 
features. Blending flux vector- and difference-splittings in the AUSM, Wada and Liou [51] 
have achieved a promising scheme, termed AUSMDV, which also cures difficulties found 
in each constituent scheme. 

In this paper, we propose another approach for upwind splitting : namely the field 
by field Hybrid Upwind Splitting. The formalism blends in a natural way the comple- 
men properties exhibited by FVS and FDS methods. Indeed, the procedure we propose 
yields upwind methods built to share the very robustness of FVS schemes in the capture of 
nonlinear waves and the accuracy of some FDS schemes in the resolution of linear waves. 
These new schemes are derived following a general hybridization technique directly per- 
formed at the basic level of the field by field decomposition involved in FDS methods. In 
particular, we stress that our hybrid upwind splitting is free of tunable parameters. 

The format of this paper is as follows. In Section 2, we introduce a convenient formal- 
ism for upwind methods that provides us with a “unified” framework for upwinding. Our 
framework turns out to be sufficiently large to recover the FDS methods and at the same 
time the FVS ones. Furthermore, its formalism enables a natural distinction between these 
two approaches. Roughly speaking, we shall see in what sense FVS methods completely 
ignore the concept of field by field decomposition which is precisely the root of all the FDS 
methods. This rather obvious distinction is revisited here for its deep consequences, lead- 
ing to in particular a very natural guideline for the design of our Hybrid Upwind Splitting 
methods. 
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Derived in full details in Section 3, these new methods fall into the framework we 
have developped in Section 2. Some of their features will be analyzed there. Section 4 is 
especially devoted to the applications to the Euler equations. Different HUS methods of 
interest are given an explicit and self-contained form of expression. Finally, we present in 
Section 5 several numerical experiments directly devoted to measuring the relevance of the 
HUS methods proposed hereinafter. 


I. Notation 


In this paper, we consider numerical approximations of weak solutions to the following 
Cauchy problem for nonlinear hyperbolic systems of conservation laws 


du dfju) 
dt dx 
u( 0, x) = uo(x). 


t > 0, x € R, 


( 1 . 1 ) 

( 1 . 2 ) 


Here, u(t, x) belongs to a phase space U € Ti p and f denotes a smooth flux function, 
f : U — ► 7 Z p , whose Jacobian matrix A(u) = df(ti) is assumed to be diagonalizable with 
p real eigenvalues Ajt(u), 1 < k < p. As is well known, this problem in general does 
not admit smooth solutions so that weak solutions in the sense of distributions must 
be considered. Moreover, an entropy-satisfying condition must be added to select the 
physically relevant (or entropy) weak solutions. We refer to Lax [28] and to Godlewsky 
and Raviart [16] for precise statements and detailed results. 

We shall assume in the sequel that (1.1) is strictly hyperbolic, that is the Jacobian 
matrix A(u) admits p distinct eigenvalues. The right eigenvectors, corresponding to eigen- 
values arranged in increasing order, are denoted by 


ri(tx), r 2 (u), ..., r p (u), 


(1.3) 


while 

l I ( u )» fJOX (l- 4 ) 

refer to the associated left eigenvectors. 

The mapping u — ► r*(u) is called the Arth field. Throughout this paper, we shall 

assume that all the fields are either genuinely nonlinear in the sense of Leix [28], i.e. (after 
a suitable normalization) 

VA*(u) • rjk(u) = 1 , Vu € U , (1.5a) 
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(1.56) 


or linearly degenerate (see [28]), i.e. 

VAfc(t/) • rjt(u) = 0, Vu € U. 

The following definitions, although standard, will be of constant use in the sequel. If the 
fcth field is gen uin ely nonlinear, a k shock wave moving at speed a is defined to be a 
discontinuity of u subject to the following Rankine-Hugoniot condition 

<t(ur - Ui) = (f(ufl) - f(«L))- (!-6) 

This k shock is said, after Lax [28] , to be admissible if 

A*(ul) — a > 0 > Ajt(ttjj) — a. (1-7) 

Associated with the Jfcth genuinely nonlinear field, a k rarefaction wave is defined to be a 
solution of the following ordinary differential equation in the phase space U 

^=r t (u(s)), s€R. (1-8) 

In view of requirement (1.5a), note that the parameter s in (1.8) is in fact restricted to 
an open subset of 71 where the mapping s -> X k (u(s)) strictly increases along the integral 
curve of (1.8) oriented by r^. 

As for the linearly degenerate field, a k contact discontinuity is defined to be a dis- 
continuity that satisfies (1.6) as well but with a — A*(uz,) = A *(«r), or a fcth wave along 
which the differential equation (1.8) is trivially satisfied. 

In the following, we shall describe numerical approximations to weak solutions of (1.1). 
Let h and h x be respectively the time and space steps. We consider piecewise constant 
approximate solutions u h : R+ x R —* U € R p ", i.e., for t € [nh, (n + l)h)[, n € JV, 

u h (t,x) = u", X e [(i - -)h x ,(i + -)h x [, i € Z, (1.9) 


l r(i+\) h * 

= — / u 0 (x)dx, i e Z. 

h x J(i-^)h x 


( 1 . 10 ) 


For n > 1, (u™ +1 ),ez is defined by the following 3-point explicit scheme in conservation 


n + 1 - (f n , - f n 


u,' = u 


(/r + x-/r-A 


h x yji + 


(i.ii) 


Here, the numerical flux defined by 



f( u ? » u T+i) 


( 1 . 12 ) 
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is assumed to satisfy the following smoothness and consistency conditions : 


/ : U x U — » TV is a locally Lipschitz — continuous function (1.13) 

/(u,u) = f(u), Vu€W. (1.14) 

By the well known theorem of Lax and Wendroff (see [20]), if the limits of approximate 
solutions (1.9), as h — » 0, exist a.e. in the sense of bounded L) oc convergence and provided 
that the assumptions (1.13)-(1.14) are satisfied, then these limits are weak solutions of 
(1.1). We refer to [16] for an exhaustive presentation and the updated references therein. 


II. A General Setting for Upwind Splitting Methods 

In this section, we introduce a convenient formalism for upwind methods. This formal- 
ism enables the classical FDS and FVS methods to be re-expressed in the same framework 
and furthermore it encompasses a new and general concept for upwind biasing (see Section 
III). 

In the sequel, we shall focus our attention on upwind schemes in the sense of the 
framework we now introduce. This framework is based on the fruitful notion of paths 
introduced by Dal Maso, Le Floch and Murat [9] to deal with nonconservative products 
arising in a hyperbolic system in nonconservation form. Following their definition, we 
consider a fixed family of paths $ in U , that is a map $ : [0, 1] x U x U — * U € W assumed 
to satisfy the following properties of consistency and smoothness, 

$(0;ul,ur) = ul and, $(1 -,ul,ur) = ur, for any ul and ur in U, (2.1 a) 
*(« ; it, u) = u, for any u in U and any s in [0, 1], (2.16) 

and 


For every bounded set K in U, there exists a constant C > 0 such that 
for any states (ul,ur) and ( vl,vr ) in K : 


d$(s;u L ,UR) d$(s;t>z,,^ft)| 


< C |(itj? — ui) — ( vr — vl)|i a.e. s in [0, l].(2.1c) 


ds ds 

As an immediate consequence of (2.1c), we have for every ul and ur in K 

d$(s\u L ,u R ) 


ds 


\< C \ ur — ui |, a.e s in [0, 1]. 


(2. Id) 


6 



We shall need also in the sequel to consider matrix-valued functions B"^ and B ■ U x 
[0,1] — + Mat('R p ,'R p ) such that : 

B + (u,s) (respectively B~(u,s)) is diagonalizable with only 

real and nonnegative A + (u,s) (resp. real and nonpositive A + (u,s) / ) 

eigenvalues for all u in H and any given s in [0,1]; (2.2 a) 


and 

For every bounded set K in U, there exists a constant C > 0 such that, 
for all u and v in K, a.e. s in [0,1], 

||B + (u, s)|| < C, ||B-(u,s)|| < C; (2-26) 

||B+(u,s) -B + (i;,s)|| < C\u - v|, ||B - (u, s) — B“(v,s)|| < C\u — v|. (2.2c) 

Based on a family of paths $ and two matrix-valued functions B ± , we propose the following 
framework for upwind biasing. 


Definition and Proposition 1 

Let $ be a fixed family of paths in the sense of (2.1) and let B + and B be two 
matrix-valued functions satisfying (2.2). $ and B ± are assumed to satisfy the following 
compatibility condition 

£ (b+($(s; u L ,u R ),s) + £“(*(«; u L , u R ), s)) d * {s ^ UR lds = f \u R ) - f(u L ). (2.3) 
Under the CFL-like condition 

t~ ,95? l A ?( u ’ 5 ^ - b 

fix. 


let us now consider the following two averaging 


2 h f 1 N ^d<f>(s;u L ,u R ) j 

ul{ul,u r ) = u l - — J B ($(s;u L ,u*),s) ^ ds, 

2 h t x \ s d$(s-,u L ,u R ) 

u R {u l ,u r ) = u r - — j ^ B + {$(s;u l ,u r ),s) ^ ds. 


(2.5 a) 
(2.5 b) 


We define a (^jB*) scheme as follows : 



? u i) + «i(«r » 


u 


*+i 


>)■ 


( 2 . 6 ) 
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Assertion . This definition yields a consistent S-point scheme in conservation form for 
which the associated numerical flux is written : 


f(u L ,u R ) = ^(f(u L ) + f(un)- 


J | B + ($(s-,u l ,ur),s) - B ($(s;mi„uk),s)} 


ds 


ds 


). (2.7) 


Rem arks 

2.1- The averaging (2.5a) may be understood as the integral conservation law applied 
over the rectangle (—h x / 2,0) x (0, h ) to the approximate solution of a specified Riemann 
problem (see the CFL-like condition (2.4)). Conversely, the averaging (2.5b) has the same 
physical picture but using the rectangle (0, h x j 2) x (0, h ). As a consequence, formula (2.6) 
can be seen as the averaging under the CFL condition (2.4) of two neighboring and non- 
interacting half approximate-Riemann solutions. Such a formal interpretation obviously 
applies to the FDS methods (we refer the reader to Remark 2.6 concerning the FVS 
schemes). 

2.2- Requirement (2.3) is an essential compatibility condition to be met by the family 
of paths $ and the B ± functions with respect to the problem (1.1). Indeed, summing 
(2.5a) and (2.5b) gives : 

^(ul + ur) = + - y J {B + +B~ }($(s; u L ,u R ), s ) d$(s, ul, ur )^ (2.8 a) 

This relation can be thought of as applying the integral conservation law over the rectangle 
(— h x /2, h x / 2) x (0, h) to the approximate solution quoted in the above Remark 2.1. Under 
the CFL condition (2.4), the compatibility requirement (2.3) clearly enforces the above 
averaging to coincide with the relevant one provided by the Godunov exact Riemann 
solver : 

^(ul + ur) = ^(ul + ur ) - ^-(f (u R ) - f(u L )). (2.86) 

In view of (2.8), the compatibility condition (2.3) might be also read as a consistency re- 
quirement with respect to (1.1). In fact, (2.3) yields a stronger requirement for consistency 
than the usual condition stated in (1.14). In particular, (2.3) is responsible for the upwind 
nature of the methods under consideration as stated in (2.9) below. Section III will give 
examples of numerical flux functions under the form (2.7) that are Lipschitz continuous 
and consistent with f in the sense of (1.14) but such that (2.3) is not satisfied. 
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2.3- For a given 5 in [0,1], neither B + {.,s) nor B (. ,s ) are assumed to be Jacobian 
matrices. Therefore, the vector-valued integrals in (2.5) and (2.7) are in general path 
dependent. 

According to D efin ition 1, when all the signal velocities associated with the numerical 
flux f{u L ,u R ) are nonnegative, i.e. B~ identically vanishes along the path $(.;ul,ur); 
then from the compatibility condition (2.3), we easily deduce that the numerical flux- 
function (2.7) must reduce to : 

f(uL,UR.) = f(uL)- ( 2 ’ 9 °) 

On the other hand, when all the signal velocities associated with the numerical flux 
f(uL,u R ) are nonpositive; i.e. B + identically vanishes along the path <&(.;u L ,u R ) then 

f(u L ,u R ) = f(u R ). ( 2 - 9fe ) 


As a consequence, we state 


Corollary 1 

Any given (<&; B ± ) scheme is a 8-point upstream method (1.9)-(1.12) in the sense of 
(2.9). 

The definition we propose meets the usual notion of upwind bias. Conversely, we do 
not claim that any given consistent and conservative 3-point upstream method in the rather 
vague sense of (2.9) belongs to the ($;S ± ) formulation. However, the examples given 
below show that our formalism is quite natural and indeed sufficiently general to recover 
the classical approximate Riemann solvers (such as flux difference splitting methods) as 
well as the flux vector splitting methods. Moreover, this formalism enables a precise 
distinction to be made between these two approaches. Roughly speaking, we shall see 
that FVS methods completely ignore the concept of waves (or field by field) decomposition 
which is by contrast the starting point of all the FDS methods. 


Proof of the Assertion 

Let us first introduce the following two functions fi and f R :Ux.U—> 7l p , setting for 
any ul and u R in U 


f 1 . .d$(s;uL,u R ) , 

fL(uL,UR) = f(«L)+ / B ($( s ; ul , ur ), s ) ^ ds, 

J 0 

f 1 ,, . . .d^(s;uL,u R ) 

f r( u Li u r) = f(ua) - J o B + ($(s;u l ,ur),s) ^ ds 


(2.10a) 

(2.106) 
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Equipped with (2.10a) and (2.10b), the averaging (2.5a) and (2.5b) can be rewritten under 
the form 

2 h 

ul ( ul , ur ) = ul ~ ^-(/l(ul,ur) — f(«L)), (2.11a) 

2 h 

ur(ul,ur) = UR - — (f (ur) - /r(«l,«r))- (2.116) 

In view of the compatibility condition (2.3), we have 

Il(ul,Ur) = /r(ul, ur) = - (/l(«L, Ur) + Ujj)) 

= /(«l, ur), Vu l ,u r eu, 


so that it makes sense to use (2.7) as a numerical flux-function and owing to (2.11), 
definition (2.6) clearly yields a 3-point scheme in conservation form. 

We now turn to proving the local Lipschitz continuity of the flux-function (2.7). Given 
a bounded set K of U, then for any (ul,ur) and ( vl,vr ) in if we write 


f(vL,v R ) - f(u L ,u R ) = ^(f(vtf) - f (ur)) + hf(v L ) - f (u L )) 




d$(s;uL,vfl) 

ds 


\B\($(s-,u l ,ur),s) 


d$(s;uL,t/fl) 

ds 


}ds. (2.12) 


Note that the above integral can be split into three terms according to 


fl + I 2 + 1$ — + g j |-^l(^(' S i u Io l; ft)5 S ){ 


d$(s-,VL,VR) d$(s;uL,VR)' 


ds 


■ds 


ds 

+ \f {|^l(^(s;VL,t?fi),s) - |F|(<&(s;ut,UH),s)|^fe^ ,t,R ^ 
1 N j5$(s;«l,«r) d$(s-,u L 

“2 1 * 


ds 


|ds. 


Using successively the smoothness condition (2.1c) satisfied by $ and (2.2b), we obtain 
the following crude estimates 


\h\<C (||£ + || + || B ||) {|(uii - v L ) - ( v R - «i,)|} 
< C \vl - ul\. 


and similarly, 


\h\ < C |u/i — ur\. 
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Moreover, in view of (2.1d) and assumption (2.2c), we deduce the following bound for J 2 

|/ 2 I < C max {| $(s;v L ,v R ) - $(s;u l ,ur) |} |vr - «lI, 

*€[ 0 , 1 ] 

< C (|ul| + \vr\) max (| $(s;vl,vr) - ${s-,ul,ur) |}. 

- Vl *€[0,1] 

Note that for s € [0,1] : 

I $(s;vl,vr) - $(s;ui,ur) | <| (*(s;v l ,vr) - v L ) - ($(s;ui,,uk) - at) | + | v L -u L \ 

Jo ^ 

< s c { I (ur — vl) ~ {ur — vl) |}+ | VL — UL | 

< Cjl V R — UR I + \ VL — UL |}? 

where we have used (2.1a) and (2.1c). Therefore, there exists a constant C depending on 
K such that 

I f{v L ,VR) - f(u L ,U R ) I < C{\ VR ~ UR | + I VL - UL I}. 

This gives the required result. 

We now conclude by proving the consistency of the numerical flux (2.7) with the exact 
flux-function f . Using the consistency condition (2.1b) concerning the path $, we have for 
any u in U 

^■ u ’l l= 0 , a.e. s in [0, 1], (213) 

ds 

since (2.1b) does not depend on s, so that /(u, u) = f(u). This concludes the proof. □ 

We now show how the classical FDS and FVS methods can be re-expressed using the 
formalism of Definition 1 . 

II. 1 Flux Difference Splitting (FDS) Methods 

Flux difference splitting methods, also referred to as approximate Riemann solvers, 
achieve upstream biasing on the basis of a given field by field decomposition. Such a 
decomposition intends to restore, but via a much simpler structure, some of the important 
features of the Godunov exact Riemann solver. Essential among the features required to 

be satisfied is the property of conservation (2.8b). 

A full hierarchy of approximate Riemann solvers exists, ranging from the simplest 
one with only two waves and one intermediate state [22] to the most sophisticated ones 
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involving a complete set of p waves separating (p — 1) intermediate states (for instance 
Roe [45], Collela-Glaz [4]). The Osher- Solomon method [35] is not strictly speaking an' 
approximate Riemann solver but nevertheless clearly falls by construction into our setting. 
We show below how to re-express any given approximate Riemann solver in terms of the 
formalism we have introduced in Definition 1. For the sake of clarity, the derivation we 
propose is performed at a formal level although it can be rigorously justified using the 
notion of graph completion introduced by DalMaso, LeFloch and Murat [9]. For this 
construction to be valid, we ask the physical space U to be a convex set. Moreover, 
approximate Riemann solutions are assumed to be made of only simple waves as for the 
exact solution (cf. assumptions (1.5a) and (1.5b)). These assumptions indeed simplify 
the derivation. Applications of this formal construction to classical approximate Riemann 
solvers are postponed to Appendix 1. 

By definition, a FDS method provides us with a self similar function w(.; ul, ur) : f € 
It — *• w(£-,ul,,ur) 6 U where f = xft, such that w(— oo-,ul,ur) = ul, w(+oo-,ul,ur) = 
ur. Generally speaking, such a function, referred as to an approximate Riemann solution, 
exists provided that ul and ur are close enough. We shall assume that the approximate 
Riemann solutions under consideration are piecewise smooth except on a finite set of 
points (£d)i<d<d m *x where u;(.; ul,ur) admits left and right limits u>(£±; ul, ur). This 
assumption is satisfied by the known FDS methods. Therefore, the total variation of such 
solutions is bounded, i.e. 

rV.j-“(u>(.;<ii,u R )) = £°° | ^( S ;«t,u R ) | d V < C, (2.14) 

where the above gradient is understood as a vector- valued bounded measure acting on It. 
For convenience, we shall adopt the strict inequality in the CFL condition (2.4), that is 
to consider times t, t < h during which the waves are strictly contained in the rectangle 
] — h x /2,h x /2[x[0,h[. It is thus equivalent to restricting our attention to f €] — fo,+fo[ 
where 

6 = I- (215) 

With this restriction, we clearly have 


TV* (w(.; u L , ur)) = TV+^ # (w(.; u L , u R )) 

since for all £ € 1Z + ; £ > fo , we have w(—£-,ul,ur) = ul and w (+£; ul , ur ) = ur. 

For an approximate Riemann solution to be relevant, we seek w(.; ul,ur) to satisfy 
the averaging consistent with the exact Godunov Riemann solver (cf. (2.8b)), 


j~ [ w(j-,u L ,u R )dx = hu L + ur)~ ^-(f(u R ) -f(u L )). 
u x J —h x jl « 2 tl x 


(2.16a) 
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(2.166) 


In view of (2.15), this requirement reduces to 

[ w(£ - , ut, ur)(1£ = £o (ul + ur) — (f( u fl) — f( u i))* 

J-S 0 

Now, let C(w) =] - 6, +to[\{Uh<d<d m „ denote the open set of points where w(.;u L ,u R ) 
is continuous. If we notice that we have : 

r A ^rn ax 

«fc(w) d=i 

/*+£ o 

Co(«l + ur)- / w(f;«L,u/i)d£, (2.17) 

</-?o 

then the relation (2.16) can be re-expressed as follows 

r A ^ m ax 

/ X] “*,«*)) = 

yC( W ) Of rf=l 

+ (f(uft) - f(«ir))- (2.18) 


The identity (2.18) plays in fact a central role in the derivation we perform below. 
Roughly speaking; when considering (2.18), the compatibility condition (2.3) mdeed sug- 
gests that the physical states u>(f; f €] — fo,+fo[, may serve as a natural variable 

to define the states involved in the path $(.; ul,ur) and that f may serve, depending on 
its sign, to define either B + or B~ . Indeed, notice that wherever the self similar function 
u>(.; ul,ur ) is smooth and non-constant in a neighborhood of £, then f returns the velocity 
at which the state w(£-,u l ,ur) propagates. Moreover, &, d € {1, ...,d max }, does coincide 
with the velocity of the dth discontinuity separating w(£J;ul,ur) from w(^’,ul,ur). 
Our purpose in the sequel is to develop a precise construction to the above assertion in 
particular regarding the treatment of the discontinuities. A geometrical picture depicted 
in Fig. 1 is useful to understand the following mathematical construction, in which Fig. 
la is a set of physical states displaying various wave structures. 

To facilitate our construction of the framework, it is convenient to introduce the 

following notation: 

TVi (o (w) = J \-^(y;uL,u R )\dy, Vf e] — £o, +fo[- (2.19) 

Let us consider the mapping X : ( € H -+ X(f) € [0,1] given by 


x(0 = 


TV+l°(w) 


TVic 0 (w(.]u L ,u R )). 


( 2 . 20 ) 
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Note that (2.20) defines a left continuous increasing function that is discontinous across 
each discontinuity shown in Fig. lb. To make use of this discontinuous data for our 
following purpose, we define a pair of functions 0 ± : £ € [— £o> +£o[ -4 ^ ± (0 €] — £o> +£o[ 

0 + (O = min {//; X(fi)>X(()}, (2.21a) 

m€]-€o.+6>[ 

and 

0~(O = min {n; X(/i) = X(0}- (2.216) 




Figure 1. Phase diagram: (a) Wave distribution, (b) Total variation, (c) 9 ± ((), (d) 
X (/x), (e) F(s), and (f) Path $. 
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Note that £ — ► |0 + (£) — 0 _ (£)| yields a left continuous step function, as depicted in Fig. 
lc. Indeed, wherever the approximate Riemann solution w(.;ul,«r) is smooth and non- 
constant in a neigborhood of ft then |0 + (ft - 0”(OI is «sro since 0 + (£) = 0~(() = ft 
Otherwise, this function gives the lenght | Ik | of the interval I* for a given k such that 

£ € Ik and 

w{h',ul,ur) = u k, VfiZlk- ( 2 . 22 ) 

Here tt* denotes the k th constant state of u>(.; ul,ur). For such values of £ € h, we indeed 
have 

9~(£) = min ft ^ + (() = max(. (2.23) 

For (2.22) to be valid, it is implicitly assumed that if there exist two adjacent intervals 
I k and Jjt+i, i.e if the constant states u k and u k +i are separated by a discontinuity, then 
the associated lenghts | I* | and | Ik + 1 I are distinct. In what follows, this assumption 
is supposed to be met. (If not, we would have to modify | 0 + - 0~ | using a relevant 
^-dependent factor to achieve a discontinuity from I k to I k + 1-) 

Let us denote by ft the location of the I th discontinuity of the function £ €] - ft, ft H 
|0+ _ 0 _ |(£). For completness, we introduce ft max +i = +ft i n suc ^ a wa y t * iat t ^ ie ^ 
U'-[^( I+1 [ yields a partition of [-ft,+ft[- In the sequel, J, denotes the interval 
[ft, ft+i [• For our purpose, we now want to remove from this partition any interval Ji such 
that | 9 + — 0 | is zero within Ji+i- 

In this way, we construct an interval S of [— £o,+ft[ using the procedure 


So 


J 0 , if |0 + — 0 - | is not zero within J\ 
0, otherwise; 


(2.24a) 


Si — { 


Ju 


if |0 + — | is not zero within J 2 

otherwise. 


(2.246) 


Note that So and Si cannot be simultaneously empty sets. Let us denote by no the minimal 
value of £ within So U Si- We then complete our construction for k < l max as follows 


k—1 

Sk = [Mfc,M*+i[, Mfc = Mo + yi 1^1, M*+i = M*+ I Sk | 

1=0 


(2.24c) 


where 


{ M, 


\s k \={ 


lo, 


if |0 + — 0~\ is not zero within Jjt+i 
otherwise. 


(2.24J) 
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Note that when | 5* 0, i.e. Sk ^ 0, then 

(~(k = Hk, V/z £ Sk, V£ € Jjt; (2.25) 

so that S* can be deduced from Jk by a suitable translation. Finally, we define 

5 = u[, moi S'fc. (2.26) 


We now reconstruct a new function X : [i £ 5 — *• X(/z) £ [0,1] using the following 
procedure. Let us consider a given nonempty subset Sk € 5, then for /z € S'* and its 
related f € <7* given by (2.25), we define 


f*(0. if0+(O-^-(O = o 

X(/z) = ^ (2.27a) 

{ X(9~(0) + n(fi)(X(9+(0) - X(d~(0), otherwise. 

Here, we have set for /z £ Sk 

?*(/*)= wher e f-£* = /z-/z*. (2.27b) 


This step does nothing but connecting the jumps of w(.; u l, ur ) with a straight line, as seen 
in Fig. Id. As a consequence, X : /z € S — ► X(fi) £ [0, 1] is a continuous strictly increasing 
function. It therefore makes sense to consider its inverse function Y : s € [0, 1] —*■ Y(s) £ S 
such that for a given s £ [0, 1], F(s) = /z where ^z is the unique solution in 5 of s = X(/x). 
Y is depicted in Fig. le. 

Equipped with this function F, we finally construct a path $ as follows for s £ [0, 1]: 


J w(£;ul,u*), if0 + (O~0 (0 = 0 

$(s-,u L ;u R ) = < , 

(O'i u L,ur) + H(O(w(0+(O'> u l,ur)-w(0 (^);?zl,ur)J; 

(2.28a) 

where for s £ [0, 1], £ is the solution : 


3 Jb€{l, Imax } ; 


( ( - 6 = ^ - fik 


( Y(s) = /z. 


(2.286) 


Notice that according to the above construction, a constant intermediate state u k separated 
by two waves from the decomposition under consideration reduces to the point Uk € 
$(.,«i, ur) while two states itk and Uk+i separated by a discontinuity (either a shock or 
a contact) are connected along $(.;«£, u.r) by a straight line (see Fig. If). Indeed for 
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s k 0 such that | 9 + — 9~ | is not zero within Jjt, there exists an index d € {1, ••-,d max } 
such that 


' w(9-(();u l ,ur ) = w(£ d -,u L ,u R ) 

<w($ + (0’,ul,ur) = w($;u l ,ur), f - Zk = V - Vk, ( 2 - 29 ) 

, * + (0 = & 

where is the location of the <fth discontinuity of w{.\ul,ur). 

The case of the rarefaction waves, that is of the states w(£; u l ,ur ) with £ €] - Co, fo[ 
such that 9 + (£) = 9~(£) = £, is straightforward. 

Now turning to the matrices B ± , we set for s € [0, 1] 


B + (<f>(s;u L ,u R ),s) = ^(0 + (O+ I ^ + (0 I )-Id-Rr, (2.30a) 

B-($(s-,u l ,ur),s) = ^(0 + (O“ \ e+ (0 I) • Id ™ , i (2.306) 


where £ is the unique solution of the problem stated in (2.28b). Note that the two matrix- 
valued functions B ± defined in (2.30) trivially satisfy the requirement (2.3) for upwind 

biasing. 

Furthermore, in view of (2.18) satisfied by the approximate solution w{.;ul,ur), we 
clearly deduce from (2.29) and definitions (2.28) and (2.30) the compatibility expression: 


/V 

Jo 


+ B- 


x ,3 $(s;ul,ur) _ 

■)($(s;ul,ur),s) -q s ds - 


Sk\ 


(■&+ 1 A w 

y, I (& u L,u R )d£ 

\9+-9- 1=0 ^ ( ' S 


dmax 

+ y U(tv(&;u L ,u R ) - w(ZJ-,u l ,ur)) = f{u R )-f(u L ). 
d=l 


From the construction performed above, the fixed family of paths $ associated with a given 
approximate Riemann solver is in fact made of the physical states involved in the wave 
decomposition in the zone of smoothness of the approximate solution w(.;u L ,u R ) while 
the matrix-valued functions B ± return the associated wave velocities. Finally, a straight 
line is used to connect in the phase space U two states separated by a discontinuity and 
B ± are defined thanks to the associated discontinuity velocity. We emphasize that the 
family of paths $ and its associated B ± functions completely depend on the approximate 
Riemann solution under consideration, that is on the given FDS method itself. 


In Section III, we shall need to make a partition of the path $(.;ul,ur) (2.28) as- 
sociated with a given FDS method, into l m ax nonoverlapping subpaths $*(.;u l ,ur) with 
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k e {l,...,/mar}- The kth subpath $ k (.;u L ,u R ) : s € [s*,s*+i] — ► $ fc (s;ti£, u R ) € U is 
associated with the kth. subinterval S k = [/zjt, of S in the following way. Here s k and 
sjt+i in [0, 1] are chosen to be 

$(•**; ^l, u R ) = w((k',UL,UR); $(sk+i;uL,u R ) == w((k+i',UL,u R ). (2.31) 

We then define 

$ fc (s; u L , u R ) = $(s; u L , u R ), s € [s*, s*+i], (2.32) 

to rewrite 

$(.;«i,u R ) = ufci-# t (.;u L ,uji). (2.33) 

For convenience in the notation, we shall make use in the sequel of 

$k(s;uL, u R ) = $jt(sfr + s(sfc + i — Sfc);ui,,ujj), Vs €[0,1], (2.34a) 

so that now 

$*(0;«l,ur) = «>(£*; u£,,ufl), $kO-)UL,u R ) = w(£ k +i;u L ,u R ). (2.346) 


Besides the formal similarity of their expression, the revisited approximate Riemann 
solvers are by construction path dependent. Put in other words, their associated numer- 
ical flux-function entirely depends on the waves (or field by field) decomposition under 
consideration. In turn, the design of such a decomposition is not only responsible for the 
merits but also for the deficiencies of the resulting approximate Riemann solver. 

Besides robustness, one of the critical issues in the derivation is how to achieve high 
resolution of stationary discontinuities while at the same time ensuring enough nonlin- 
ear dissipation to select the physically relevant discontinuous solutions. The difficulty in 
achieving this goal comes from the fact that entropy violation only occurs for stationary 
or near stationary discontinuities associated with the genuinely nonlinear fields. 

The current prominent approximate Riemann solvers have been designed to perfectly 
resolve a stationary contact discontinuity ([4], [35], [45]) and to yield a sharp steady discrete 
shock [4] or to exactly capture a shock ([35], [45]). But associated with this desirable 
accuracy, their nature of approximation results in insufficient robustness in calculating 
some practical problems. This is in particular true for the Roe and Collela-Glaz methods 
that may fail in the capture of large expansion waves, both admitting nonphysical shocks 
as stable solutions. A lack of robustness has been reported for the Osher scheme in the 
capture of strong detached shocks [24]. 
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Fixes have been introduced to cure some of the discussed failures [12], [21], [24], [46], 
but at the expense of losing accuracy. These extra mechanisms need to be properly timed 
up either to be efficient or to restore locally some desirable accuracy that was lost by the 
uncorrected scheme. This is, in particular, the case of the widely used entropy fix proposed 
by Harten [21] which is applied in practice to all the fields of the Roe decomposition in 
order to increase the overall robustness of the method. This procedure in turn destroys the 
high resolution of the Roe scheme for stationary discontinuities. We stress that the exact 
capture of a contact discontinuity at rest is a crucial requirement for the correct resolution 

of boundary layers [50], [14] (see also section V). 

We now turn to re-express FVS methods in the framework of Definition 1. 

II.2 Flux Vector Splitting (FVS) Methods. 

In this section, we study a different and widely used approach for implementing upwind 
biasing, namely the FVS methods. The statements given below will underline the deep 
differences from the FDS methods we have just revisited in our setting. We first state 

Definition 2 

Let (/ + ,/ - ) be a pair of continuous and piecewise differentiable functions f . U 
Tl p , required to satisfy the following two conditions 

f + (u) + f-(u) = f(u), ueM, (2.35) 

and in any region of smoothness 

B + (u) = df + (u) (respectively B~{u ) = df~(u )) is diagonalizable 

with real and nonnegative (resp. real and nonpositive) eigenvalues. (2.36) 

For any u L and u R m U, the FVS method defines the numerical flux-function associated 
with the pair (/ + ,/~~) to be 

f(uL,u R ) = f + (u L ) + f~(u R ). (2.37) 


Rem arks 

2.4- The exact flux-function f needs not be a homogeneous function of u of degree one 
for a consistent FVS method to exist (in the sense of (2.35)), as pointed out in [22]. 

2.5- In view of the requirement (2.36), the definition of the FVS methods we give is 
different from the one originally given in [47] and [22]. However, we stress that all the 
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FVS methods that have been brought to our attention, satisfy the condition (2.36) (see 
Lerat [29] for the Steger- Warming splitting [47]). Moreover, the Van Leer method satisfies 
by construction the above definition while by contrast it does not fall into the framework 
for splitting methods proposed in [47] (see Hirsh [23]). The implications we give of the 
Definition 2 will be assessed in what follows (see, in particular, Propositions 2). 

As an immediate consequence of Definition 2, we state 
Corollary 2 

Let be given a general FVS method associated with the pair (/ + ,/“). This FVS 
method is a (4>; B ± ) method for any given fixed family of path 4*. Here, B ± are the 
Jacobian matrices defined by (2. 86). The numerical flux function defined in (2.7) can be 
re-expressed as 

/(«!,««) = \ (f(ut) + f(«R) - J\b+ - B-)(»(«;m.uit)) a * (s; <k), (2-38) 

or under the equivalent form 

f(u L ,u R )=^(f(u L ) + f(u R )- Y sign(±l)(f ± (uR)-f ± (uL))). (2.39) 

2 -i,+i 


FVS methods therefore fall into our framework for upwind biasing but according to 
our formalism, FVS methods are path independent. We stress that by this strong property, 
a general FVS method is in complete opposition with the basically path dependent FDS 
methods. This deep distinction, as provided by our formalism, will be given a physical 
interpretation in the next section devoted to the Euler equations. 

Rem ark 

2.6- A formal explanation of the path independence can be given on the basis of the 
following two systems of conservation laws : 


du df (u) = 
dt dx 

du <?/+( «) 

dt dx 


associated with the following initial data 


u(0, x) 


k ur(x), 


if x < 0, 
if x > 0. 
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We shall ass tim e, although it is not clear with no further assumption on / + and / , that 
both the Riemann problems above admit solutions, repectively denoted by u>(., u L , u R ; f ) 
and w(.,u L ,u R ;f + ). Now, in view of Remark 2.1, the averaging (2.5a) applied to 

w{.,UL,UR\f~) yields at least formally 

2/i 

u L (u L ,u R ) = U L - ^-(/ _ ( ur )-/"(« l )), 

since under the requirement (2.36), no wave can propagate outside of the rectangle 
(~h x / 2,0) x (0, h). Conversely and using again (2.36), we can write for the averaging 

(2.5b) but with u>(., ul, ur; f + ) 

2/i 

Ur(ul,UR ) = UR — — (/ + (mr) — / + ( U l))- 

Notice that the above averaging can be rewritten under the expected forms 

ul(ul,ur ) -u L - |“ ((/“(“«) + / + («l)) - f («L))i 
Ur(ul,Ur) — UR — j—(f(uR) — (/ (ur) + / + (“£-)))’ 

where we have used the consistency condition (2.35). 

The averaging procedures performed above have made no use of the precise wave 
patterns involved in w(.; u L , ur] /*). The averaged solutions would have been the same 
for any given other wave decompositions, provided that the waves under consideration are 
only propagating in their respective relevant rectangle (see requirement (2.36)). The path 
independence property stated in Corollary 2, may be understood in that way. 

To conclude the present remark, we note that FVS methods can be given at least 
formally a Godunov- type interpretation (but in sharp contrast with FDS methods), under 
which they make use of two distinct sets of Riemann problem for evaluating the same 
numerical flux (cf. in particular Remark 2.8). 

Proof of Corollary 2 

The path independence property clearly comes from the Jacobian nature of the ma- 
trices B ± . Indeed, we have for any given fixed family of paths in U 

f — <*» = (2 ' 40a) 

J B-(*(»;uL,KR)) °*^’^ [ t ’ t * — <fa = (2.406) 
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Therefore, substituting (2.40a) and (2.40b) in (2.38) clearly gives the expected flux formula 
(2.39). Now using the consistency condition (2.35), let us write 

f(u L ) + f(u R ) = (f+iiiL) + f~(u L )) + (f + (u R ) + f~(u R )). (2-41) 

Formula (2.37) directly follows by substitution in (2.39). 

Moreover, by differentiating the consistency relation (2.35), we deduce the following 
identity 

B + (u) + B~(u ) = A(u), a.e. u € W, (2.42) 

where A is the Jacobian matrix associated with the flux function f. So the compatibility 
condition (2.3) stated in Definition 1 is satisfied for any given path 4>. This concludes the 
proof.D 
Rem arks 

2.7- Although the identity B + (u ) + B~(u ) = A(u) holds true for almost every u in 
U , we do not have in general B + (u) — A + (u) and (therefore) B~(u) = A - (u) for all u 
in U. (see Steger and Warming [47]). If these equalities were met for all u, the associated 
FVS method would be necessarily path dependent since A ± (u) are not Jacobian matrices, 
except precisely when all the eigenvalues of A(u) keep the same sign. Notice that even for 
such u, the equality B ± (u) = A ± (u) is not necessarily satisfied. This is, in particular, the 
case of the T-Boltzmann schemes developed in Appendix 2. 

2.8- Generally speaking, B + (u ) and B~(u) do not commute. As a consequence, they 
do not admit the same set of right eigenvectors. In that case, notice that the right eigen- 
vectors {rjfc(u)}i<fr< p of A(u) do not diagonalize either B + (u) or B~(u). This property is 
indeed responsible for some difficulties in the analysis of general FVS methods. 

By contrast with the property stated here, formulae (2.30a) and (2.30b) clearly in- 
dicate that the B ± matrices do commute for a given FDS method; furthermore, FDS 
schemes satisfy an averaged consistency condition (2.16) while FVS methods merely obey 
the poiniwise consistency condition (2.42). To our knowledge, only the Osher scheme satis- 
fies a similar condition (note that Van Leer [49] has recognized the Osher scheme as a FVS 
method in the case of a single conservation law; for a closely related result, see Brennier 
[ 2 ]) 

These observations further underline the fundamental differences between FVS and 
FDS methods. 

By the property of independence with respect to path, we have 

Proposition 2 

A general FVS method in the sense of Definition 2 cannot exactly capture a stationary 
discontinuity. 
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Proof of Proposition 2 

Let us consider two states u l and ur connected by a single discontinuity at rest, either 
a shock or a contact discontinuity. In view of the averaging (2.5a) and (2.5b), this wave is 
exactly captured if and only if 

ul(ul,ur) = ul and ur(ul,ur) = ur, (2.43) 


which is equivalent to asking for 

[ = (2.44) 

for all path $ connecting ur to ur. Taking the scalar product of (2.44) with ^ 
and integrating over t € [0, 1], we necessarily have 


f 1 f 1 d$(s;uL,UR) , d $( s ; ur,ur) 

dt - d ^~ L B + ($(s-,u l ,u r )) ^ ds+ 

jr - ,u P )\ di{s '^' U ^ dsdt = 0 . ( 2.4b) 


By assumption, for all u in U, B + (u ) only admits nonnegative eigenvalues. So that, we 
necessarily have 


ds os 


(2.46) 


for any given path $. Taking advantage of the path independence property, we can always 
choose a path $ such that 


9$(s; Ut ,U R) B+Ws; Ui UR)) 9f(£^K ) > 0i s e I C [0, 1), (2.47) 

ds us 


where I denotes a nonnegligeable open subset of [0, 1]. Otherwise, B + would be identically 
zero for almost every u in U\ this would in turn imply, thanks to (2.42), that B — A does 
admit strictly positive eigenvalues (since A does), which contradicts with our assumption 
(2.36). So the strict equality (2.47) holds true for some paths 4>. Notice that from a generic 
smoothness assumption (2.1d) made on for any u L and u R in U we have 


a $( 3 ; ML ,« g )| = 0(\ur-u l \), a.e. s €. [ 0 , 1 ], 
os 


Therefore, we choose paths $ such that 

[' u L , UR ))^^ds > - u L \\ 


(2.48) 


(2.49) 
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where S is a positive constant independent of ul and ur. On the other hand, writing for 
a smooth path, for a.e. s € [0, 1] : 


d $( t ; uL , UR ) 5$(3;u£,,ur), c^S^sju^ur) , „ 

—^1 1 < max — (ttR-ttL \\s 

at as s€[o,i] a z s 




{A.OU) 


Notice that the normalization has been accounted for in the last parenthesis. We see that 
if the path under consideration is chosen with the following additional smoothness 


d 2 $(s;u L ,UR) 

d 2 s 


= 0(|ur - uz,|). 


a.e. s € [0, 1], 


(2.51) 


we necessarily have 



d$(t-,u L ,u R ) 

dt 


d$(s-,UL,UR). + d$(s-,UL, Ur) 

95 + ($( s ; m l , u r )) dsdt = 


0(\u R — tt£,| 3 ). (2.52) 


So that the identity (2.44) is clearly impossible if ul and ur are chosen close enough 
while connected by a single discontinuity at rest (notice that this is always possible). This 
concludes the proof.Q 

Proposition 2 applies to states u i and ur connected by a single shock wave with zero 
speed or a contact discontinuity at rest as well. This result tells us that by contrast with 
the FDS methods revisited in section 2.1, a stationary contact discontinuity is necessarily 
smeared when approximated by a general FVS method. In fact, in the applications, it 
turns out that such a linear wave is unendly smeared, in contrast with a nonlinear wave as 
discussed above. The numerical dissipation underlined here has an important consequence 
in practical issues when deeding with the numerical approximations to some viscous form of 
problem (1.1) and (1.2). Its implications have been in particular emphasized by Van Leer 
et al. [50] for the numerical approximation of the Navier-Stokes equations, thus making a 
general FVS method an irrelevant candidate to this purpose. 

The proof given above makes the expression (2.45) play a central role. This expression 
can be related to the amount of the numerical dissipation inevitably performed by a FVS 
method (for the relations between the rate of the Entropy dissipation and the numerical 
viscosity we refer to the work by Coquel and Le Floch [5], [6]). Its quadratic strength 
is in complete agreement with the results obtained in [5] for numerical methods whose 
numerical viscosity never vanishes (see [6] for the deeply different results obtained in the 
case of FDS-type Godunov method). According to this observation, the Jacobian nature of 
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B ± combined with the relevant assumption (2.36) on the signs of their eigenvalues therefore 
seems to inevitably yield numerical methods with nonvanishing numerical viscosity. 

Put in other words, any attempt to design pure FVS methods with vanishing viscosity 
or even with very low diffusivity should result in a violation of the stability requirements 
(2.36); namely B ± would necessarily admit eigenvalues with the wrong sign. The issue is 
a marginal stability for the method and at least a loss of monotonicity (at first order of 
accuracy) for such low diffusive FVS methods. These conclusions are in complete agree- 
ment with the observations by Liou and Steffen [31] and Coirier and Van Leer [3] and make 
“pure FVS-method a dead-end street” for viscous calculations [3]. 

Besides this negative result, the robustness achieved by FVS methods strongly suggests 
that some of their advantageous features should be kept in the derivation of new upwind 
biased procedures. With respect to the formalism provided by Definition 1, FVS methods 
are entirely characterized by their associated matrix-valued functions B ± as it is asserted 
by Corollary 2. As a consequence, B ± provide us with the essence to be kept from the 

FVS method. 

Moreover in order to get high resolution of stationary contact discontinuities while 
preserving robustness, Proposition 2 clearly suggests that at least extra information specif- 
ically regarding the contact discontinuity must be recognized. Since the failure met by the 
pure FVS methods inherits from their path-independence property, the hybridization we 
have to perform must result in a path-dependent method. Therefore, this additional and 
essential step must relate (to a large extent) to a fixed family of paths. As revealed by 
the fr am ework we have introduced earlier, natural families of paths already exist in the 

current FDS methods and they can be readily utilized. 

All these observations naturally suggest that the tools to be involved are a combina- 
tion, within the framework of Definition 1, of the two existing approaches for upstreaming. 
The following section describes a way for hybridizing these tools in order to achieve the 
requirements that have been put forward so far. 


III. Field by Field Hybrid Upwind Splitting (HUS) Methods 

The concept of hybrid upwinding we propose, finds a natural definition in the setting of 
($ ; B±) schemes we have introduced in the last section. This new approach for upwinding 
is achieved by combining the main interesting features of two distinctly different concepts. 
The motivation for such an approach stems from several practical and theoretical reasons 
outlined in the Introduction. Specifically, we want to derive a new class of upstream 
methods with the following properties : 
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( PI)- keep the very robustness of flux vector splitting in the capture of genuinely 
nonlinear waves (strong shock as well as large rarefaction waves), 

( P2)- preserve the property met by the FVS methods to select only entropy satisfying 
approximate solutions, 
while at the same time 

( P3)- offer the decisive ability of some flux difference splitting methods to exactly 
capture stationary contact discontinuities. 

In view of requirement (P3), we shall only consider in the sequel ($; Bp DS ) approxi- 
mate Riemann solvers that perfectly resolve a contact discontinuity at rest. 

As suggested above, HUS methods are derived by making a given FVS scheme path 
dependent with respect to a family of paths coming from a given FDS method. To this 
end, we notice that the first two requirements (PI), (P2) only address genuinely nonlinear 
fields, while by contrast, requirement (-P3) is only related to linearly degenerate fields. 
As a consequence, for a fixed family of paths $ provided by a given FDS method, we 
are naturally led to describe each “genuinely nonlinear” subpath with FVS- type Bp VS 
Jacobian matrices while keeping unchanged for “linearly degenerate” subpaths the original 
Bpos ma trices that are compatible with $. 

With clear notation, GNL(<&) denotes the part of $ made up of the subpaths associ- 
ated with the genuinely nonlinear fields while LD($) stands for the remaining subpaths. 
For our following purpose, we restrict our attention to considering FDS methods that 
satisfy the following requirement LD($): 

E f ( b fds + B FDs)(^k(s; u L , u R ), s ) ^2f^iLL^l ds = 

1 fc; ♦*€ LD($)' / ° 

E (f(^Jk(l; UL, UR )) - f($*(0; UL, Ur))) . (3.1) 

k ; $*€ LD($) 

This we shall refer to in the sequel as a compatibility condition for hybridization. Remark 
3.1 will assess such a reference. The validity of (3.1) is investigated below for some of the 
prominent FDS methods; examples and counter examples are given. We show in remark 3.2 
that requirement (3.1) indeed provides us with an essential condition for ensuring upwind 
biasing for the hybrid methods we now introduce. 

Definition and Proposition 3 Let be given a general (. ,Bp VS ) FVS method and a 
($, Bp DS ) FDS method subject to the condition (8.1). A Hybrid Upwind Splitting (HUS) 
method is defined by the following consistent numerical flux-function : 

f(uL,U R ) = ~ (f(ui,) + f(uj?) 
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^ Z 1 ,, , ^d$k(s;u L ,u R ) 

- 5 Z / |B|FVs($Jt(s;WL,«R» ds 

$ fc €GATL($r 0 

$ fc € LD($r° 

The numerical flux (3.2) can be equivalently revrritten under the form 

f(u L ,u R ) = f F vs(uL,u R ) + | X) / (I B If^s - |B| F Ds) a$fc(5 ^ L,U — ( 3 * 3a ) 




or 


/(uL,UR) = /FDs(ttL,ttR)-| E / (l^lfV5-|B|FD5) a$fc(S ^ L — *■ ( 3 ‘ 36 ) 


<fr*eGiVL($) 


Assertion A HUS method is a ($; B±) scheme in the sense of Definition 1. 

The relevance of our definition for a hybridization technique will be assessed in what 
follows with regards to theoretical as well as practical issues. In particular all the HUS nu- 
merical flux-functions, derived below in our new upwinding setting, preserve the Lipschitz- 
continuity of the underlying FDS and FVS methods. For simplicity in the notation, the 
endpoints of a given subpath (.;u L ,u R ) are from now on denoted by 

$*(0;ul,ur) = = «*+i- ( 3 * 4 ) 


Proof of the Assertion 

We have to check that the hybridization procedure introduced above yields upwind 
methods in the sense of Definition 1. We therefore have to show that the compatibility 
requirement (2.3) satisfied by both the underlying FDS and FVS methods is preserved 
when exchanging for a given subpath, i.e., the B ± matrices, according to the nonlinear 

nature of the associated field. 

In this way, let us consider 


r 1 d$(s]UL,u R ) , 

I (Bhgs + B hus )($(s-, u l,ur),s) s ~ 

V 0 

is . ULyUR)) ?^^Sl ds+ 

Jo 


y: f (Bfvs + B FV s)(^ k 

$ k €GNL(*) J ° 


(3.5) 


4 > k €LD (< !>) 
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The first set of integrals related with the GNL subpaths can be rewritten as 

(Bfvs + — ds - 

E A B fvs + b fvs) ($* («; u L , u fi )) 1 ^ > ' Ur) ds, (3.6) 

where in view of the pointwise consistency property (2.42) met by the FVS methods, we 
have 

J ( b fvs + Bp VS )(^s]u L ,u R ),s) ^^j-^^ ds = f(«/t) - f(u L ), (3.7) 

and 

E / ( B FVS + B FVs)(^k(s]U L ,U R ))^^^^ds = 

<f> k €LD(*) J ° 

E (f(tifc+i)-f(«*))- (3-8) 

4> k €LD(*) 

Now from the identities (3.6) to (3.8), taking into account our assumption (3.1) for the 
underlying FDS method clearly enforces the relation (3.5) to reduce to 

J o ( b hus + B hus)(*( s ’’ u “*)* ~ ds = f ( u «) “ 
which is the required result 

Rem arks 

3.1- As clearly shown by the identity (3.8), the condition (3.1) comes from the use of 
FVS-type B ± matrices along the GNL subpaths of the FDS-path $. Such a condition is 
clearly required for the compatibility condition (2.3) to hold true for the HUS methods. 
Corollary 1 thus applies and ensures upwind biasing in our framework for hybridization. 
(3.1) therefore turns out to be a compatibility condition for a given FDS method to provide 
a suitable candidate for field by field hybridization with FVS schemes. 

3.2- In order to stress the relevance of the compatibility condition (3.1), let us assume 
here that it is not satisfied by the underlying FDS method. We show below that an attempt 
for hybridization with any given FVS scheme according to (3.2) yields a method that fails 
in upwind biasing. Let / denote the resulting numerical flux. We introduce the following 
formal averaging under the CFL condition (2.4) 

2 h 

ul{u l ,ur) = ul- —(f{u L ,u R )-f(u L )), (3.9a) 

2/i 

ur(ul,ur ) = UR — — (f (ur) — f(v,L,UR)). (3.96) 
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Simple algebraic manipulations in (2.5a) and (3.9a) yield: 

/* 1 'U'Lj ^ ft} 

2 (f{uL,UR.) — f( u L)) = % Yj / BpYg^k^ULj'UR)) ^ 

/^irr /A\ ^ 0 


GiVL($) ' 


+ 2 2_^ I B FDS (^k( s 't u Li u R.)-> s ) Q s 

LD($r° 

+ E /’ (fi?vs + B; vs )(» t (.; “1.. »«)) 9 * t (i a“ — ^ 

E f ( b *fds + (310) 


LD(*) ' 


Thanks to the identity (3.8), we see that the formal averaging (3.9a) does not coincide 
with the relevant one stated in (2.5a) unless the requirement (3.1) is satisfied; it clearly 
exhibits unexpected contributions from waves supposed to travel outside of the domain. 
To go further, ass um e now that for the states u l and ur under consideration, LDi^ $ ) 7^ ® 
and that all the states connecting ul to ur in IA via 4>(.; ul, ur) satisfy 

B fvs ($(s-,ul,ur)) =0, Vs € [0,1], (3.11a) 

B fds ($(s;ul,ur),s) = 0, Vs G [0,1]. (3.116) 

That is both the FVS and the FDS methods under consideration give for such u L and u R 


/fVs(«L,«h) = fFDs(uL,UR.) = f(«L)‘ 


(3.12) 


Then, in view of (3.11), formula (3.10) clearly yields 


2 (f(u L ,u R ) - f (u L )) 


Y (f(«k+i) - f ( u *)) 

LD($) 


/) / B FDS (^k{s',ULy u R)T S ) q s aS 

■ . 0 


LD{$) 


7^0. 


Without the condition (3.1), f(u L ,u R ) does not return the expected flux (3.12). Such a 
result holds true for any given pair of states ul and ur that satisfy (3.11) with LD { $) 7 ^ 0* 
It is easy to transpose this negative conclusion to (3.9b) to see that conversely when all 
the signal speeds are nonnegative, / does not return in general the expected backward flux 
f(uft). Such features result in a failure for upwinding. As a consequence, we may expect 
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the method, which is not an HUS method according to our definition, to suffer from a 
lack of stability. Changing the underlying FVS method is obviously of no help. These 
developments enlighten the essential role played by the compatibility condition (3.1) for 
hybridization. 

3.3- Since the requirement (2.3) holds true for FDS methods, the condition (3.1) is 
equivalent to asking for 




(f(tlfc+i)-f(«*))- 

$*€GNL($) 


(3.13) 


Indeed, note that we have 

^2 (f(ujt+i) -f(«*» + ^2 (f(uk+i)-f{uk)) = f(u R )-f(u L ). 

$t€GAri,($) €LD($) 

Before investigating for validity the compatibility condition (3.1) for some of the 
prominent FDS methods, we first state two straightforward consequences of the Defini- 
tion 3 : 

Corollary 3 

If ul and ur are connected in the phase space with only linearly degenerate subpaths, 
that is GNL($) = 0; then any given HUS numerical flux-function (S.2),(3.8b) reduces to 

f(uL,UR) = fFDs(uL,UR). (3-14) 

For a single contact discontinuity at rest, we have in particular 

f{u L ,u R ) = f{u L ) = f(u R ). (3.15) 

Thus, requirement (P3) is met. Conversely, we also clearly have 

Corollary 4 

If u L and u R are connected in the phase space with only genuinely nonlinear subpaths, 
that is L£>($) = 0; then any given HUS numerical flux-function (3. 2), (8. 8a) reduces to 

/(«l,«r) = fFVs(uL,u R ). (3.16) 
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In view of the equivalent flux formula (3.3a), the second term on the right hand side may- 
be understood as an antidiffusive procedure to correct a general FVS method with respect 
to requirement (P3) as stated in Corollary 3. This antidiffusive flux is designed to exactly 
capture a stationary contact discontinuity and, roughly speaking according to Corollary 4, 
to vanish for genuinely nonlinear waves in order to restore the ability of the FVS methods 
for capturing large expansion waves and entropy-satisfying strong shocks. 

Note that in the same way, the equivalent flux-formula (3.3b) may be conversely 
understood with an additional diffusive procedure designed to enhance the robustness of 
the underlying FDS method in the capturing of nonlinear waves and/or to enforce the 
selection of physically relevant approximate solutions. Section V devoted to the numerical 
experiments will in particular assess this important issue. 

The above interpretation in terms of an additional antidiffusive flux when considering 
formula (3.3a) or a diffusive flux, but with respect to (3.3b), must be given a more precise 
foundation. For such an interpretation to be valid, we should at least expect the following 
inequalities to be satisfied for ui and ur in U sufficiently close : 


[' £ [\\ B | F VS - I B | F D s) ^ t( 7/’ UR) <fa} > 0. 

f U L , U R ))dt{ Y, [ (\ B \FVS-\B | FDS) ^ 1 “ 0 ’ 

Jo to [ 'GNL(*) Jo 


where U uu denotes the Hessian of given entropy U associated with (1.1). The investigation 
of these two requirements is however beyond the scope of the present paper and will be 
addressed in a forthcoming work. 


For the current FVS and FDS methods, numerical evidences (see Section V) clearly 
demonstrate the relevance of the hybridization procedure we propose and strongly support 
the fact that both the linear and the nonlinear stability requirements (P1),(P2) are met 

as well. 

We now turn to studying the validity of the compatibility requirement (3.1) for some 
of the prominent FDS methods : namely the Osher- Solomon, the Collela-Glaz and the 
Roe method. They are derived as a (^,5 ± ) scheme and will be revisited in Appendix 1 
under our formalism. We first state 


Lemma 1 

Both the Osher- Solomon and the Collela-Glaz methods meet the compatibility condition 
(2.1) for hybridization. 


31 



Proof of the Lemma 1 

For some of the basic features in use below, we refer the reader to Appendix 1. We 
consider at first the Osher-Solomon method. By construction, for any given subpath 
€ LD($), there exists <xjt € H such that the endpoints ujt = $*(0 ',ul,ur) and u*+i = 
<I>jfc(l; tt£, ur) satisfy the following Rankine-Hugoniot relation 

<nk(ttjfc+i - u k ) = (f(«* +1 ) - f(tt*)). (3.17) 

Here, cr k is nothing else but the velocity at which the fcth contact discontinuity travels. 
Now since the method returns (see Appendix 1) 

J (B+ + = °*( u H-i ~ u k ), (3.18) 

the identity (3.17) clearly ensures the requirement (3.1) to hold true for the Osher-Solomon 
method. 

Again by construction, the Collela-Glaz method does satisfy both the equalities (3.17) 
and (3.18) for a given <r k € 72,. The condition (3.1) is thus also valid for this scheme. □ 

As a consequence, both the Osher-Solomon and Collela-Glaz methods provide us with 
relevant candidates for a field by field hybridization. By contrast, we show below that a 
Roe solver does not necessarily meet the requirement (3.1) for a suitable hybridization. 
Indeed when considering a Roe type linearization B{ur , u^), standard arguments yield for 
a given € LD($) (see again Appendix 1) : 

f ( B + + «-)(*»(.; u L ,u R ), s) ^^Sl ds = AiOii.uaXufcH - »»). (3-19) 

where A t(ut, ur) denotes the associated eigenvalue coming from the Roe type linearization. 
Here and by contrast with (3.17), (\ k (uL,v,R);u k ,u k +i) is solution of the following jump 
relations 

A k (uL,UR)(u k+ i - u k ) = B(uL,UR.)(u k +i — Ufc). (3.20) 

Note that we only have a priori 

T Afc(«L, UR)(u k+1 - u k ) = f (u R ) - f(tt L ) - y A fc(uL, u R )(u k+ x - u k ) 

&k£LD($) $*€GWL($) 

since B{ul,ur)(ur — ujf) = f(uj?) — f(u/,). With no further assumption on the Roe 
type linearization B(ur, ur), there is no reason for the above expression to reduce to the 
expected identity (3.1). Indeed, a failure for satisfying (3.1) can be reported for the original 
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Roe method devoted to the Euler equations. This failure is inherent to the fact that the 
associated intermediate states do not share in general the same pressure and the same 
velocity. Thus we have checked that 

Lemma 2 A general Roe type linearization does not necessarily yield a FDS method 
that satisfies the compatibility condition (S.l) for hybridization. 


Rem arks 

3.4- In view of the identity (3.17), we see that both the Osher-Solomon and the Collela- 
Glaz methods actually obey a stricter compatibility condition than the one required m 
(3.1). Indeed notice that we have by contrast with (3.1) 


r 1 , _l w _ , \ ,d$ fcQs;ui„UR) , _ 

I {&FDS + -^FDs)(^ fe ( S > UL ' 1 UR )‘> 5 ) Q s 

Jo 

>)(M 5 


/ 


( Bfvs + Bfvs) 


= f(ut+l) ” {< - Uk) ' (s21> 


for any given € LD($). 

3.5- In a recent work devoted to the Euler equations, Mehlman [33] has interpreted the 
Collela-Glaz method in terms of a Roe solver. Such an interpretation makes nonempty the 
set of Roe-type linearizations for satisfying the requirement (3.1). Up to our knowledge, 
the Collela-Glaz method seems to yield the only yet known Roe method that satisfies (3.1) 
but in a stricter way owing to (3.21). This distinction might leave room for the existence 
of Roe type linearizations that only obey the weaker form (3.1). 


In the sequel, we shall assume that there exist Roe type linearizations with the fol- 
lowing additional consistency condition 

V B(u L ,un)(u t+1 -u t )= £ (f(u i+ i)-f(«)), (3.22) 

* k €LD { $) *fc€LD(*) 

which is nothing but the compatibility condition (3.1) for hybridization. Such an assump- 
tion makes the derivation of Roe-type hybridizations a rather formal issue. Nevertheless, 
such HUS methods exhibit some discrepancies of practical interest with the ones coming 
from both the Osher-Solomon and the Collela-Glaz schemes. Assumption (3.22) is actually 
made to provide one with a better insight into the hybridization procedure we propose and 
to highlight some of its main features. 

We now turn to giving a closed-form expression to the numerical flux-function (3.1) 
resulting from a hybridization of a general FVS method with the FDS methods we have in- 
vestigated : namely, the Osher-Solomon or the Collela-Glaz methods and the Roe schemes 
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but subject to (3.22). The statements below indicate the simplicity for formulating a 
numerical flux for all of these hybridizations of interest. 

Concerning the Osher-Solomon and the Collela-Glaz type hybridizations, we have 

Proposition 3 

Let be given a general FVS method defined by the pair (/ + ,/“). For any given ul 
and ur inU, the numerical flux function associated with the ($os; Bpvs,os) HUS method 
(respectively the ($CG',Bp DSCG ) HUS method) reads 

/(«I,XR) = / + («t) + /■(«*)+ Y, «(/“(«*+!) -/“(“*))> (3.23a) 

k ; $*€££>($) 

with 

e k = -sign(<Tk). (3.236) 

Here, a k = Aj k (u k ) = A*(u*+i) denotes the velocity associated with the kth linearly degen- 
erate subpath of the Osher-Solomon path (resp. the Collela-Glaz path) (cf. (3.17)). 

In a similar way; assuming the existence of relevant Roe methods for hybridization, 
we state 

Proposition 4 

Let be given a general FVS method defined by the pair (/ + ,/ - ) and a Roe-type lin- 
earization B(ul,ur) subject to the additional consistency condition (3.22). For any given 
ul and ur in U, the numerical flux function associated with the ($Roe‘, -^rvs, Roe) HUS 
method admits the following form 

f(u L ,UR) = f + (u L ) + f~(u R )+ €jfc|(/ € *(ujk +1 )-/ £ *(ufc)) 

*; **€LD($) '' 

+ i^A*(uL,Uft)ujt +1 - f(u fc+ i)) -(X k (uL,u R )u k -f(tt*)) j, (3.24a) 

with 

e k = — sign(A k (u L , ur)) . (3.246) 

Here, \ic(ul,ur) denotes the eigenvalue of B(ul, ur) associated with the kth “ linearly 
degenerate ” subpath of the Roe path (cf. (3.20)). 

Remarks 

3.6- The Osher-Solomon and Collela-Glaz type hybridizations look formally the same. 
Only, the endpoints ujt are different. 

3.7- By contrast, the Roe- type hybridization admits a somewhat different form of 
expression. However, this expression is deeply related to the previous ones if we observe 
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that the Rankine-Hugoniot 
(3.17) 


relations for a k contact discontinuity ensure when considering 


CTjfcUfc — f(ujfc) = CTjfcUfc+l — f(u*+i)- 


Of course, this equality does not hold true for the Roe method which on the other hand 
satisfies the following distinct property of continuity 


- B(uL,UR)Uk = Ajt(«I,MR)«Hl “ B(v.L,UR)Uk+ 1- 

already stated in (3.20). These equalities are in general incompatible, except precisely 
when the kth. contact discontinuity under consideration is the only wave separating 
and ur. 

3.8- Note that despite the additional requirement for consistency (3.22) to be met by 
the Roe methods, we do not have in general 

y; ejtj (Afc(uL,UR)ujt+i — f(ufc + i)) — (Afc(uL,u,R)ufc — f(u*:))} =0, 
fc ; i t SLD ($) 

except precisely when 

ejt = ej, VAr, /; $*, € LD($), 

that is when all the eigenvalues associated with the linearly denerate subpaths share the 
same sign, as it is expected in the light of Remark 3.2. 

3.9- In view of the above remarks, the discrepancies seen in the formulation of the 
HUS numerical flux-functions come from the way by which the compatibility condition 

(3.1) is satisfied by the underlying FDS method. According to the remark 3.4, both the 
Osher- Solomon and the Collela-Glaz schemes obey the stricter condition (3.21) while the 
Roe type linearizations under consideration are only expected to satisfy the weaker form 

(3.1) . 

Let us stress the fairly simple form resulting from the numerical flux-functions ob- 
tained so far. By comparison with the flux formulae (A1.5), (A1.10) and (A1.16) given in 
Appendix 1 , these expressions do not involve the endpoints of any of the inner genuinely 
nonlinear subpaths, as opposed to all the original approximate Riemann solvers. We also 
emphasize that sonic points for rarefaction waves play no role in our Osher-type hybridiza- 
tion. As a consequence, the hybrid numerical flux-functions proposed herein are really 
simpler than their underlying approximate Riemann solvers. 

This important practical issue will be further exemplified in the next section devoted 
to the applications to the Euler equations. Indeed, we shall appeal to the symmetry 
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property of f ± with respect to the Mach number to directly evaluate f €k and not to deduce 
f~(u) from / + (u) through the consistency relation (2.35) as it is commonly done. These 
observations allow another practical interesting aspect of our new approach for upwinding. 


Proof of Proposition S and f. 

We introduce the following function <j> = 2(fnus — fFVS ) 7Z P . In view of 

the flux formula (3.3a), the antidiffusive flux <j> reads 

<f>(uL,UR)= ^2 [ (I^Ifvs — |-B|fds) — -- ds. (3.25) 

LD(*) J ° L ° 

Let us first consider the Osher-Solomon type hybridization. We can rewrite <f> with clear 
notation as: 

4>(ul,ur)= ^2 <t>k(uL,UR ), (3.26a) 

k ; €££>(*) 

where 

4 > k { uL , UR ) = j (|#|fvs — |B|f£>s ) — ~ — d $. (3.266) 

Taking advantage of the Jacobian nature of the Bp VS matrices, we easily get 

/ |B|Fvs(^it(-s;uL,UR)) ds - 

(/ + (« *+i) ~ /+(«*)) - {f-(u k+1 ) - f-(u k )). (3.27) 

On the other hand, thanks to Appendix 1, we have 

J \B\FDs($k(s; ml, ur, _| Gk | ( Ufc+1 _ Uk ) 

= sign(oie){f(u k+ i) - f(u k )}. (3.28) 
Using the consistency relation (2.35) for the FVS method under consideration, we rewrite 
f (u k+ i) - f(u k ) = {f + (u k+1 ) - f + (u k )} + {/-(ujt+x) - /"(a*)}, (3.29) 

so that (3.27), (3.28) and (3.29) yield 

' -2 (/“(ujfc + i)-/~(ujk)), if a k > 0, 

4>k{uL,UR) = < (3.30) 

k +2 (/ + (it fc+1 ) - f + (u k )), otherwise, 

which is the required result. 
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It is easy to see that this derivation, also applies to the Collela Glaz hybridization 
along exactly the same lines. 

We now come to the Roe-type hybridization. Following the lines of the derivation 
above, we rewrite <j> as 


^(ul,«r)= <f>k (ul,ur), 

k; * k €LD(*) 


(3.31a) 


with 

<f>k(uL,UR) = [{f + ( u k+l) ~ f + ( u k)} ~ {f ( u k+l) ~ f ( u *)}) 

- \^k{uL,UR.)\(uk+i - Uk )• (3.316) 


Here, we have used the fact that 

J \ B{uj j ,Ur) IfDS = l B{ u Li u r) |fD5 ( u k + 1 — Uk) 

=| A r(«l,«r) I ( u k + 1 — Uk). (3.32) 

To make the function <f>k given by (3.31) similar to the one obtained in (3.30), let us write 

I A*(ul,«k) | (u k +i - u k ) = sign(A k )({A k (u L ,UR)iik+i - f(u k+ i)} 

— {^k(uL,un)uk — f(ufc)} + {f + {uk+i) ~ f + ( u k)}+{f (uk+i) — f («*)})> (3.33) 

so that (3.31), (3.32) and (3.33) give the expected result. This concludes the proof.D 

We turn now to study the Lipschitz-continuity of the HUS methods we have considered 
so far. 

Proposition 5 

Let be given a general FVS method defined by the pair (/ + ,/ )- For the Osher- 
Solomon, the Collela- Glaz and the Roe type hybridizations, the associated numerical flux - 
function respectively given by (3.23) and (3.24) is a locally Lipschitz continuous function 

uxu-*n p . 

Proof of Proposition 5 . 

Let us again consider the antidiffusive flux <f> = 2(fnus — fFVs) :U xU —> Tl p given 
by (3.3a). Note that fHUS and 4> have the same smoothness. 

First, in the case of the Osher- Solomon type hybridization, we recall that <f> can be 
written as 

4>{ul,ur)= *^2 <j>k(u L ,u R ) (3.34a) 

fc; $*€LDO) 


37 



with 

( -2(/“(u*+i) - /“(«*))» *7 > 0, 

^*(«L»“ft)=\ • (3.346) 

{ +2(/ + (u*+i) - / + (u*)), otherwise. 

In view of the underlying ODE (1.8) (see also (Al.l)), the endpoints Uk(uL,v,R) define 
smooth mappings IA x H —*■ U. As a consequence, we only have to study the continuity of 
<f>k in the open set D* = {u l ,ur € U\ |cr*.| < e} with e arbitrary small. By substracting 
the two equations in (3.34b), we get 

2(/ + (u* +1 ) - /+(«*)) + 2 (/-(«*+!) - f~(u k )) = 2(f (u* +1 ) - f(«0). 

Here we have used the consistency relation (2.35) : f(«jt) = f + ( u k) + f~( u k)- For a 
stationary contact discontinuity, i.e. <r* = 0, we observe that f(ttfc+i) — f(txfc) = 0. This 
proves the required result. 

It is easy to see that this conclusion also holds true for the Coilela-Glaz type hybridiza- 
tions. We recall that for this approximate Riemann solver, the endpoints ujt(tii, ur) define 
locally Lipshitz continuous functions. 

We now come to study the smoothness of the Roe-type hybridization. Following the 
same lines as above, we have to study the smoothness of <f>k(uL,UR) given, for positive 
values of Ajt(ui,, ur), by 

<f>k = — 2(/ - (u*+i) - /“(«*)) 

- ((A(itL, UR)uk+i - f(ujt+i)) - (A h(ul, UR)uk - f(ufc)) , (3.35a) 

and for negative values of A jt(ttL, ur ), by 
<t>k = 2(/ + (u fc+ i) - f + (u k )) 

+ ((\(uL,UR)uk+i - f(ujt + i)) - (Ajt(«L, Uft)u* - f(«*)) • (3.35 b) 

We shall assume that the endpoints u k(u r , u r) associated with the Roe method under 
consideration, yield locally Lipshitz functions. Therefore, we are again led to study the 
smoothness of this type of hybridization for A k{uL,UR) in the neighbourhood of 0. Sub- 
stracting again both equalities (3.35a) and (3.35b) yields 

2(f(u*+i) — f(u*)) + 2^(A(ui,,t/ft)ujfc+i - f(ujt+i)) 

- (Ajfc(«L,«fi)“fc - f(«Jfc)) = 2\(u L ,UR)(uk+i - u*). 
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This ensures the local Lipshitz continuity of the Roe-type hybridization when 
A* (u£,,ur) = 0. This concludes the proofO 

Rem arks 

3.10- In the proof given above, the enforced consistency condition (3.22) to be met 
by the Roe-type linearizations plays no role. As a consequence, numerical flux-functions 
of the form (3.2) but performed with any given Roe method yield consistent f and lo- 
cally Lipschitz continuous functions provided that the Roe method in use meets this last 
property. However, we again emphasize that without (3.22) and (3.1) the resulting HUS 
method is not upwind biased. 

3.11- In the case of the Roe- type hybridization, it is essential to use the original 
Roe type linearization to describe the linearly degenerate fields. Exchanging formally 
the endpoints of the the flux formula (3.23) with the ones coming from the Roe path $, 
would yield a slightly simpler numerical flux (cf. expression (3.24)). However, we stress 

that such a function is discontinuous , since in general we have f(u*+i) - f O*) # 0 when 
\ k (u L i ur) = 0 (cf. Remark 3.6), except when u L and u R are connected in the phase space 
through a single contact discontinuity at rest. Despite this last property, this simpler 
procedure (which is not a hybrid method) must not be considered. We underline that the 
lack of continuity makes not valid the conclusion of the Lax-Wendroff theorem. 

These results of Lipshitz-continuity further assess the relevance of our setting proposed 
for hybrid upstream schemes. In fact, some refinements can be introduced in the fight 
of the previous results. Indeed, we can modify Definition 3 in order to derive simpler 
Hybrid numerical methods in the case of a path $ made of at least two distinct linearly 
degenerate fields. Let us focus our attention on the Collela-Glaz method for instance, for 
which the associated approximate Riemann solution is univalued. We propose the following 

procedure: 

a- For q, € LD($), find q 0 such that A, 0 = min, |A,| (notice that cr q = \ q (u k ) = 
A,(ufc+i)). 

b- Define the following numerical ($; B ± ) method 

For k, 1 <k <p with p ± qo, $fc( 5 ! u L , ur) is endowed with Bp VS matrices, 

and 

For k = qo,$k(s]u L ,u R ) is endowed with A . 

According to this definition, we get 

£ f \B\ F vs 94 ’ k{S '2 L — d ‘ = <■**<.**) ~ f + M) - (/■(“»> - f > l)) 

k&o C 
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+ (/ + (« R ) - / + (<x,„ + ,)) - (/-(«*) - r («»+,)). 


On the other hand, we have 


J o ds = sign(A qo )(f(u qo+1 ) - f^J). 
Standard arguments yield the following numerical fluxes 


f(uL,UR ) 


/ + 0l) + / (««)-(/ (u go+ i)-/ (u 9o )), t/A, o >0, 

, f + ( u L ) + /“(«r) + (/ + (u 9o+ i) - / + (u ?0 ))> otherwise. 


We stress that a stronger compatibility condition than that stated in (3.1) is now required 
for the above procedure to yield an upwind method : namely the condition (3.21). This 
condition is precisely satisfied by the Collela-Glaz scheme in use. Also this HUS method is 
still Lipschitz continuous and satisfies requirement (P3) while being slightly simpler than 
the one proposed above. On the other hand, this simplication obviously comes at the 
expense of an enhanced numerical viscosity. Although easy, the extension to the Osher- 
Solomon path yields a different numerical flux since the associated approximate solution is 
multi-valued, so that more than one linearly degenerate wave may travel with zero speed. 

Moreover, even simpler combinations of FVS and FDS methods can be performed if, 
in opposition with the Definition 3, the genuinely nonlinear subpaths of $ are not sys- 
tematically endowed with FVS-type B ± matrices. But it turns out that the resulting 
simplification always seems to come at the expense of a discontinuous numerical flux func- 
tion. Once again, we stress that the lack of continuity makes such method an irrelevant 
candidate for practical purposes. In the light of the remarks of this section, the definition 
we propose appears to be a simple as well as relevant setting for hybrid upstream methods. 


IV. Application to the Euler Equations 

The present section is devoted to applying the general hybridization procedure pro- 
posed in Definition 3 to the Euler equations modeling a calorically perfect gas. Detailed 
flux-formulae are given, coming from a hybridization of a general FVS method with either 
the Osher-Solomon or the Collela Glaz FDS schemes; these later schemes being revisited 
in Appendix 1 in terms of our formalism. Let us recall that such FDS methods indeed 
are relevant for our purpose since they both meet the compatibility requirement (3.1) for 
hybridization. 
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Besides the celebrated Van Leer flux vector splitting [48], we shall also focus on the 
attractive Boltzmann schemes recently introduced by Perthame [38], [39]. Based on the 
general framework for Boltzmann schemes he proposes, we derive an one-parameter family 
of such methods, deduced from a previous study by Kaniel [26]. 

The system of 2D gas dynamic conservation laws projected in an arbitrary coordinate 

x becomes 

(4.1) 


du df (u) 

1 — = 0 , 

dt dx 


where the unknown u and the flux-function f read 


u = 



f(u) = 


pv 

pv 2 +p 
pvw 

,(pE + p)v. 


(4.2) 


Here the thermodynamic closure equation is given by 

p = (7 - 1 ){pE - \p(v 2 + TO 2 )), 7 6)1,3], 7 = constant. 


(4.3) 


In the sequel, we give a closed-form expression to various HUS methods based on some 
of the prominent approximate Riemann solvers. The pair of split fluxes (/ + , / ) refers to 
a general FVS scheme. We shall restrict our attention to FVS methods that satisfy the 
following properties of symmetry with respect to the Mach number. In particular, both 
Van Leer’s and Perthame’s (see Appendix 2) split fluxes (/+,/“) possess this symmetry 
property. All other state quantities set constant, and M — c being the speed of sound, 
we ask, after Van Leer [48], for 


f}{M) = - f;(-M ), f+ (M) = f- (-M), 

= (4.4) 

Note that the exact flux-function (4.2) satisfies such a set of symmetry properties. 

IV.l HUS Methods Based on Compression/Rarefaction Waves Decomposition 


Here, we focus our attention on the Osher-Solomon approximate Riemann solver. 
For the sake of completeness, some of its basic features are recalled in Appendix 1. The 
Osher path $(.;u L ,u R ) (either in the natural or reverse ordering) connects u L to ur m 
the physical plane (t,x) using only rarefaction or compression waves. This path is known 
to exist in the large as long as there is no cavitation and it is made of respectively two 
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genuinely nonlinear and two linearly degenerate subpaths for the system (4.1)-(4.3) under 
consideration. The latter LD subpaths are in fact associated with the contact discontinuity. 

Let us denote u* L and u R the left and right states separated by the contact discontinuity 
travelling with the velocity v* . Let us recall that both the pressure and the velocity v stay 
constant across this contact discontinuity. States u* L and u* R axe determined by solving the 
following usual set of equations made of the Riemann invariants of the system considered. 


Pl P* 2 2 

7i = eV’ WL = W ' 1 ’ n ~ e (7^T) ct = v * “ e CT7) ci ' (4 ' 5a) 

Pr p* 2 2 

7r = K' wr = w *’ '' B + £ (7rTj c « = '’*+ £ (7nj c «’ ( 4 - 54 ) 

where c = ^ (4.5b), we set e = +1 (respectively e = —1) for a path in 

the natural order (resp. reverse order). The above system of equations is known to admit 

a unique and explicit solution (when no cavitation) given by the following algorithm well 
suited to our purpose: 

1 ) Compute Sl = Pl(pl ) 7 , Sr = pr(pr) 7 , and then Si = 

2 ) Compute the sound speed 



(cr + cl) 4- ^27^(t>R — vl) 

1 + Si 
2 


and deduce c R = Sic * L , 

3 ) Compute v* = v L - (c L - c * L ), 

4 ) Compute k = L c £ 2 and n = (Slk) (, ~ l) and deduce p* L = kSl, Pr = 7r Sr 

5 ) Deduce p* — np^. 

Notice that the above algorithm always returns positive pressure and densities pro- 
vided that the left and right states meet this property. Moreover besides exponentiations, 
a fairly small number of divisions axe required. 

With this solution, we apply the flux formula (3.6a) to the Euler system to get 


Lemma 1 The HUS numerical flux-function based on a rarefaction/ compression 
waves decomposition reads 


/(«z.,«r) = / + (ul) + / (ur) + 


-(/-(»%) ~ r(ol)). 

,+(/ + Mt)-/ + M))> 


if v* > 0, 
otherwise . 


(4.6) 
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Moreover , the switch involved in (4-6) can be automatically accounted for , provided that 
symmetry conditions ( 44 ) are satisfied by the underlying FVS method. In that case, the 
HUS flux reads componentwise 


U («L, «*) = /;(Ml) + f;(M R ) + (/;(- \ m r\)- rn - 1 *2 1))’ ( 4Ja ) 


f P v(u L ,u R ) = /+(Ml) + /„(Mr) - sign(v*)(f+(- | Mr |) — f+( I M£ |)), (4.76) 
fpw(uL,U R ) = Ipw(Ml) + fpw(M R ) + (fpw(— I Mr I) — fp W {- I M L I)), (4.7c) 


f pE (uL,u R ) = f+ E (M L ) + f pE (M R ) + (f+ E (- I ATr D — fpE(~ I W* (4 ' 7d) 


This significant simplification over the original Osher-Solomon scheme is brought about 
by the present hybridization procedure. Moreover, no sonic points are involved so only 
four fluxes need to be evaluated instead of six, and this combination involves no switch. 


Proof of Lemma 1 

Starting with the flux formula (3.23a), we rewrite the summation over the two linearly 
degenerate subpaths $ 2 (u£,u*), $ 3 (u*,u^), where u* denotes the intermediate endpoint. 
Note that these two subpaths axe associated with the same eigenvalue v , we therefore 
have e 2 = €3 . As a consequence, the intermediate state u* plays no role since 


E * 


(/<>(«*+,) - /“)(“*)) = «((/*(«*) - /'(”*)) + (/'(“*) - /'(<))). 

= «(/*(<*) - /'(<)). (4 - 8) 


which yields the required flux formuls. (4.6). 

Now, thanks to the symmetry properties stated in (4.4), focusing our attention on the 

first component of the flux (4.7a), we have, setting all other state quantities constant, 
if v* > 0, clearly M£ = v* f c* L =\ Mf | and the same holds true for M£; so that 

-(f~(M R ) - f;(MD) = (/+(- \Mr\)~ /+(- I Ml I)), 

This identity also holds true for the two last components : namely for pw and pE. 

Concerning the momentum equation, note that both the Mach number M L — v jc L 
and M* r = v* / c%, , share the same sign with v*. As a consequence, if v* > 0, we can write 

in view of (4.4): 

- f-(MD) = -(/*(- I Mil)- f+(- I Mi D) 
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To conclude, when v* < 0, we notice that = — | | and M R = — | M R |; therefore 

the formula (4.7b) is obviously equivalent to (4.6) in that case. 

IV.2 HUS Methods based on Shock/Shock Waves Decomposition 

Here we address the hybridization of the Collela-Glaz Riemann solver with a general 
FVS scheme. As for the Osher path, the present path is made of two genuinely nonlinear 
subpaths but associated with shock waves and of two linearly degenerate ones related to 
the contact discontinuity. As mentioned in Remark 3.5, such an hybrid method looks 
formally the same as those studied in IV. 1. Only the endpoints are different. But unlike 
the rarefaction / compression waves decomposition, no explicit formula is available and an 
iterative procedure is required to evaluate the corresponding states u* L and u* L . For the 
sake of completness, we give below an efficient algorithm for such a purpose. 

This algorithm, devised by Gottlieb and Growth [18], solves the problem in considering 
the following two functions of the velocity v : Pi(v) and Pr(v) defined as follows: 

Pi( v ) = PL + 1PL [~ — — ]Wi, 

CL 




(4.96) 

and in a symmetric way 



Pr( v ) = PR + 7 PR [- 

c>«' 

(4.10a) 

where 



w« = ^[^] + (i + { 

(4.106) 


The functions p* L and p* R respectively account for the pressure jump at a given velocity v, 
across the left and right shock waves. 

These two shock waves can connect ul and u R via a contact discontinuity provided 
that there exists a value v* of the velocity such that 

p£(»*) = Pr(»*)3P*- (4.11) 

Note that the pressure p* standing for the pressure at the contact discontinuity may be 
negative. In (4.11), v* is nothing but the velocity at which the contact discontinuity travels. 
With these notation, the algorithm therefore reads 

find v* € 71, such that Pr(v*) — Pl(v*) = 0. (4-12) 


(4.9a) 
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A Newton iterative procedure is proved to be really efficient in solving (4.12); the initial 


where we have set 

VL =Vl + 


* ZVL + VR 

V °~ 1 + Z 

. H- 1 ) 

(4.13) 

2 2 

_ CR fPL\ 

(4.14) 

(7 _l) Ct ' VR UR (7-1) R ’ 

CL 'PR' 


Once the velocity v* and hence the pressure p* are known, densities are given by the 
following relations 


Pl = PL 


(v* - Vl¥ 
(p* -pl) ’ 


Pr = PR 


(v* - vr¥ 
(p* - pr) ' 


(4.15) 


while the tangential velocities are simply w* L = w L , w* R = wr. The states u* L and u R 
are now entirely characterized, the HUS flux-function based on the shock/shock waves 
decomposition is given by expression (4.6) stated in Lemma 1. 


IV .3 On Some Boltzmann FVS methods 

In the hybrid flux formulas derived in the last sections, the pair of split fluxes (/ + , /" ) 
refers to a general FVS scheme. Besides the celebrated Van Leer FVS method, the Steger- 
Wanning splitting can be chosen also. The new and attractive Boltzmann- type schemes, 
considerably improved recently by Perthame in [38], [39], provide one with other pairs of 
split fluxes ( f + ,f~ )• 

Perthame’s motivation was to build entropy-satisfying Boltzmann schemes with fi- 
nite speed of propagation, corresponding to velocity repartition functions with bounded 
support. Besides the important entropy-satisfying property, the use of velocity reparti- 
tions with bounded support brings a dramatic improvement in accuracy over the schemes 
based on the physical Maxwellian distribution e~*~ (see for instance [10], [43]). Moreover, 
Perthame has proved under natural conditions concerning the velocity repartitions that 
both the density and the pressure remains nonnegative [38]. 

In an independent work, Kaniel [26] has designed such velocity functions with bounded 
support. These functions yield an one-parameter family of velocity distributions, formally 
related to the isentropic Euler equations written for a calorically perfect gas with P €]1, 3]. 
This family turns out to fall into the framework for Boltzmann schemes proposed by 
Perthame (see Appendix 2). The splitting is achieved on the equivalent Mach number 
yJfpTp. The key point is that any given r-Boltzamnn scheme is naturally made consistent 
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with the exact flux-function (4.2) for which in general 7 is different from I\ We only 
consider in this Section the two simplest members corresponding to T = 2 and T = 3. The 
latter case coincides with a scheme already considered by Perthame [38] and corresponds to 
the monoatomic one-dimensional gas dynamic. The choice T = 2 is related to the shallow 
water equations and to our knowledge, has not been considered by other authors. 

These FVS schemes are defined by the pair of split fluxes (/ + , /") we now describe. 
Let us set 

c 2 = — ; M = (4.16) 

P c 

where c and M respectively stand for a sound speed and a Mach number. Under the CFL 
like condition 


the pair (/ + , / ) is given by 



i\v\ + J—)<l, 


(4.17) 


if M > 1, f + (u ) = f(u), / (u) = 0; 
if M < 1, f + (u ) = 0, /“(«) = f(«); 


(4.18a) 

(4.186) 


and for | M =| f |< 1 



2 - Boltzmann 

3-Boltzmann 

c 

i/f 


|M| 

| M | 3 

|M| 2 

ft 

±|(2 ± 3 M+ | M | )pc 

±i(l±2M+ | M | )pc 

f± 

J pv 

/± V + i(3±M(4-|M|))p 

f±v+i(2±M(3-\M\))p 

J pw 

fpW 




(§t; 2 +Aj)/±±^ 


i(12 ± 45M + M 2 (40 - 7 | M |))pc 

^(1 ± 4M + M 2 (4- | M \))pcj 


(4.18c) 


A = 

The flux formula given above are expanded for the sake of similarity. Notice however 
that the 3- Boltzmann scheme can be rewritten under the equivalent but very simple form 
(see [38]) 


ft = ±f (M ± l) 2 , ft, = ± x ) 3 ' 

ftw= w f?’ 


pc 


f%.±^ M ±Xf + xlft (4.19) 
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We emphasize after Perthame that the 3 -Boltzmann scheme coincides with the Van Leer 
splitting when 7 = T = 3. Besides the simple expression, these pairs of split fluxes provide 
us with FVS methods of interest. Let us indeed recall that under the CFL condition (4.17), 
such FVS methods preserve the positivity of both the density and the pressure [38]. 


V. Numerical Experiments 

In this section, we investigate the performance of various HUS methods for accuracy 
and robustness. Intending to illustrate the various properties we have stated, we applied 
the HUS methods to several test cases ranging from stiff shock-tube problems to 2-D high- 
speed viscous flows, some of them involving strong shock-wave/boundary layer interaction 
which leads to flow separation. Such cases provide a convenient benchmark for testmg 
both the shock-capturing capabilities of a given method and its accuracy m the prediction 
of surface quantities, in particular of the heat transfert rates. Moreover, comparisons of 
the results are performed with some of the underlying pure FVS and FDS methods in 
order to enlighten the benefits of the present hybridization technique. 

Five distinct sets of problems are addressed to evaluate the relevance of the HUS 
methods. The first set deals with shock tube problems designed to be a severe test for 
shock capturing schemes. Results are assessed with the aid of an unusual but enlightening 
sensitive tool : the entropy balance equation. We point out how this additional equation 
can be used, numerically speaking, to gain a valuable insight into a general flux splitting 

method. 

The four other sets of problems are to measure the relevance of the HUS methods for 
Navier-Stokes calculations : namely the one- dimensional conical flow problem proposed 
by Van Leer, the two dimensional shock wave/boundary layer interaction experimentally 
investigated by Hakkinen et al., a Mach 25 viscous flow over a blunt body including grid 
refinement studies and finally two Mach 10 flows past compression ramps with 15° and 

20° wedges. 

V.l Results for the Shock Tube Problems 

The problems discussed here are designed for a chemically frozen mixture of calorically 
perfect gases, diatomic and monoatomic oxygen in the present set of problems; such a 
modeling makes the specific heat ratio 7 depend on the frozen mixture composition and 
therefore enables distinct but constant 7 for the left and right states. 
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Three distinct HUS methods are tested : the hybridization of the Van Leer FVS scheme 
with the Osher-Solomon FDS scheme evaluated using the exact Riemann invariants, and 
the 2- and 3-Boltzmann FVS schemes but with their associated r-Riemann invariants. 

A grid of uniform 200 points is used for all test cases. The second order of accuracy, 
when needed, is achieved using the MUSCL approach written in primitive variables. The 
CFL condition has been set at the constant value 0.5. 

Case 1 

Left state : Xoi — 1, Xo — 0, P — 9.88 x 10 5 , T — 2438, v = 0, 

Right state : Xoi = 0, X 0 = 1, P = 9.93 x 10 3 , T = 2452, v = 0. 

This first test case is designed to show that the use of improper Riemann invariants in the 
T-Boltzmann schemes does not affect their ability to correctly capture strong shock and 
rarefaction waves. Montagne et al. [34] have reported noticeable errors in the location (and 
velocity level) of the shock wave in the approximate solutions. To make the problem stiffer, 
we have deliberately chosen the initial mass fractions in such a way that jl = 1.4 and 
7 r = 1.6667 so that a small bump is observed in the velocity profiles, as shown in Fig. 2. 
The approximate solutions given by the HUS-F-Boltzmann schemes axe quite comparable 
with the others, thus confirming the relevance of the approach. A closer examination of all 
results suggests that the HUS with Van Leer and Osher-Solomon hybridization compares 
the best with the exact solution, specially in the velocity profile. 

Case 2 

Left state : P = 6.50 x 10 5 , T = 2442, v = 0, 

Right state : P = 1.00 x 10 3 , T = 346, v — 0. 

The mass fractions are calculated under the assumption of chemical equilibrium. The 
approximate solution calculated by the shock-shock Collela-Glaz approximate Riemann 
solver at the physical time t = 7 x 10 -5 s is used as the initial condition for this test 
case. This initial data presents a large nonphysical expansion shock located at x = 0.5. 
Our purpose for presenting this problem is to investigate the entropy stability of the HUS 
methods (hence only looking at the first order accuracy). To this end, we calculated the 
numerical entropy source term for a given consistent scheme according to the following 
formulae 

2/i 

ul(ul,ur) = ul- ^-(/(ul,mr) — f(ux,)), (5.1a) 

*2h 

ur{u l ,ur) = ur- — (f (u R )~ f(u L ,u R )), (5.1 b) 
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and 


^(5(ul) + S(ur) - S(ul ) - S(ur )) 

-f — (F(ufl) — F(ul)) = Entropy source term (5.2) 
hx 

for the convex entropy pair (5 = plog{ e 1 < x °*' x 2 i) t F = Su). The derivation of this 
formulae is postponed to Appendix 3. As it will be made clear there, formula (5.2) is the 
entropy equation averaged over a shifted computational cell, centered at each interface 
(CFL < 0.5). It accounts for the rate of entropy dissipation/production of a given 3- 
points scheme for a single jump (« i , ur). Figure 3 displays the velocity profile at the final 
physical time t = 2.1 x 10 -4 s. All the methods have correctly dealt with the nonphysical 
shock expansion, except the exact Roe method by Abgrall. The entropy source term is 
nonpositive, including at the glitches. The plot of the entropy disspation at x = 0.5 gives 
a measure of the expected rapid smearing of the inadmissible shock with the time itera- 
tions. Both the original and HUS Van Leer methods dissipate entropy at the same rate 
as the Osher does while the FVS and HUS T-Boltzmann schemes yield higher numerical 
viscosity — neither B + nor B~ admits a vanishing eigenvalue (no glitches are noticeable). 
We note that any additional viscosity resulting from the energy splitting due to the hy- 
bridization formulation only acts in the GNL fields, as clearly shown by the following sets 
of calculations. 


Case S 


Left state : P = 3.54 x 10 4 , T = 450, u = 0, 

Right state : P = 3.49 x 10 4 , T = 2500, v = 0. 

The mass fractions are calculated under the assumption of chemical equilibrium. The 
exact solution is determined by a slowly moving contact discontinuity at speed v* = 3 
m/s. Our goal is to assess the relevance of our hybridization technique by investigating 
its n um erical diffusion, “inside” a slowly moving contact discontinuity (if exactly at rest, 
this diffusion vanishes). To this purpose, the entropy source term defined by (5.2) provides 
a very efficient tool because its magnitude is closely related to the amount of numerical 
dissipation associated with the method. Figure 4 shows the entropy source term at the 
location of the contact discontinuity versus time, clearly demonstrating that the numerical 
dissipation (as measured by the entropy source in (5.2)) of all the HUS methods m the 
linearly degenerate fields are almost same as that of the Osher scheme and the Roe scheme 
by Abgrall. By contrast, the dissipation entropy rates for the pure FVS methods are about 
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30 times higher ! The smearing of the density profiles gives a qualitative measure of the 
amount of dissipation. 

Case 4 


Left state : P = 5.73 x 10 2 , T = 199, v = 2200 (M = 7.8), 

Right state : P = 2.23 x 10 4 , T = 546, v = 0. 

In this case, we assess the performance of numerical schemes for resolving a shock moving 
slowingly against a large pressure jump (roughly 400 to 1 ratio). Roberts [44] pointed out 
that a numerical noise could be generated at the shock by some numerical schemes, leaving 
an oscillation with frequency of several grid points. Figure 5 shows such an oscillation from 
the Roe and Osher schemes. Both HUS and Van Leer schemes give monotonic resolution of 
the trailing shock, but a slight bump at the contact discontinuity is seen to be assocaited 
with the latter. It is interesting to note that although all schemes capture the trailing 
shock well with only one to two cells, they perform rather poorly on the leading shock. 

Case 5 


Left state : P = 5.73 x 10 2 , T = 199, v = 4100 (M = 14.6), 

Right state : P = 5.73 x 10 2 , T = 199, v = -4000 (M = -14.5). 

This case provides a severe test on the capability for resolving shocks resulting from two 
opposing hypersonic flows. Figure 6 indicates that the Roe, HUS, and Van Leer schemes 
give very accurate and similar results, except there are slight differences near the initial 
contact point and at shocks. Both the natural and reverse orderings in the Osher paths 
for the HUS scheme give comparable results. 

V.2 Viscous Conical Flow Problem 

Proposed by Van Leer et al. [50], this test case provides an ideal benchmark for 
testing the accuracy, or conversely the numerical dissipation of a numerical flux-function, 
for viscous calculations. We solve at the first order of accuracy the one-dimensional conical 
Navier-Stokes equations for a Mach 7.95 viscous flow with 7 = 1.4 over a 10° cone, using 
64 cells. The Reynolds number is Re = 4.2 x 10 5 and the Prandtl number is set Pr = 0.72. 
Approximate solutions for the temperature, pressure and normal velocity are plotted in 
Fig. 7 and compared with those obtained by the Roe method. An excellent agreement 
with this benchmark scheme is achieved for all the HUS method with a slight superiority 
for the less dissipative Van Leer building block, but only in the capture of the shock wave 
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(GNL field). By contrast, the pure FVS methods predict considerably broadened boundary 
layers. These results clearly show the relevance of the hybrid method we propose for viscous 
calculations. For completness of evaluating the proposed HUS method based on the Van 
Leer splitting, we include the convergence test. The present and Roe methods are taking 
a virtually identical path to convergence, while the Van Leer method is about 50 % faster. 

V.3 Shock Wave/Boundary Layer Interaction Problem 

This test provides a comparison of the solutions obtained by the present HUS Van 
Leer method and Roe scheme with the data measured by Hakkinen et al. The two- 
dimensional thin-layer laminar Navier-Stokes equations axe solved for the flow of Mach 
2 0, Re = 2.96 X 10 5 an d Pr = 0.72. Second-order accurate solutions are carried out on a 
75 x 65 grid. As shown in Fig. 8, both HUS and Roe solutions axe comparable in accuracy. 
They give very good agreement with the surface pressure data, while the HUS result is 
seen to be slightly closer. However, both predict too big a separated region and a long 
relaxation in the reattaching region. This seems to be a common feature also reported 
by other authors; but the cause of this discrepancy is yet unknown. Unlike the conical 
flow case, we discover some interesting features in the convergence history for the HUS. 
While almost identical in the first four decades of the residual reduction, the HUS method 
departs from the Roe method by taking a steeper convergence rate, roughly a factor of 

two. 

V.4 Supersonic/Hypersonic Blunt-Body Flow Problems 

The first problem in this category deals with a high-speed laminar flow past a hy- 
perboloid configuration. The curvature radius at the stagnation point is R = 0.44 m, the 
asymptote angle is given by a = 29.3°. The flow parameters are Moo = 25, Too — 192.34 K 
and T wa ii = 800 K. The Reynolds number based on the curvature radius at the stagnation 
point is set to Re R = 7684, derived assuming a standard atmosphere at H = 77 km. The 
Prandtl number is set at Pr = 0.72 and the ratio of specific heats is given by 7 = l- 4 - 

Four grids were used for the grid-convergence study. The most refined grid, made of 
50 x 112 nodes, was generated using an exponentially stretched spacing in the streamwise 
direction while the clustering in the body normal direction is defined using the following 
ratio parameter 1 + esin( 7 r(j - l)/(nj - 1), see Gazaix [15]. From this grid, three subgrids 
50 x 56, 50 x 28 and 50 x 14 axe extracted systematically to study the influence of the 
refinement in the body normal direction. Figure 9 displays the grid plotted in every other 
interval in both directions. 
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Computations with the present HUS method and the underlying Van Leer scheme are 
compared in terms of surface quantities prediction. Calculations with the Osher method 
were also carried out on the two finest grids to provide reference predictions. Figure 10 
shows the results from the Van Leer method. The pressure coefficient (Fig. 10(a)) is not 
sensitive to grid resolution. In sharp contrast, the skin friction coefficient (Fig. 10(b)) is 
heavily affected by increased grid resolution and differs by a significant amount over the 
entire surface from a grid to another. Mesh resolution has the same drastic effect upon the 
Stanton number (Fig. 10(c)) : these surface quantities are grossly overpredicted on coarse 
meshes and significant discrepancies still exist for the two finest grids, indicating that the 
solution has not yet achieved convergence for the Van Leer method with the finest mesh. 
These results again confirm the excessive numerical diffusivity of the Van Leer method, 
already reported by other researchers [14,30,50]. This superfluity in diffusivity finds its 
roots in the lack of resolution for linear waves (such as contact discontinuity) inherent in 
the FVS schemes. 

The surface quantities with the present HUS method are plotted in Figure 11. Con- 
trary to the Van Leer scheme, no sensitivity with grid resolution is observed for all the 
surface quantities. The coarsest and finest meshes yield no significant variations in either 
the skin friction coefficient or the Stanton number. For these two grids, the hardly ob- 
served discrepancies are by several orders of magnitude less than the one existing for the 
FVS quantities computed with the two most refined grids. These results therefore indicate 
the dramatic improvement in accuracy brought by our hybridization procedure. Besides 
this significant improvement, the present method exhibits the excellent shock capturing 
capability from its underlying Van Leer scheme. Indeed, the smoothness of the predicted 
surface quantities strongly supports that none of these calculations exhibits any tendency 
towards the development of anomalous solutions, such as the carbuncle phenomena re- 
ported for the Roe method [45] or other methods [14]. 

Turning to the Osher scheme, Figures 12 and 13 depict the predicted skin friction 
coefficient and the Stanton number compared with the results from the Van Leer and the 
present HUS methods. The original Osher and HUS schemes predict virtually identical 
surface quantities for both meshes. This excellent agreement shows that the present hy- 
bridization procedure returns the accuracy of the Osher scheme in resolving the boundary 
layers. 

Another problem of interest concerns with the so-called “carbuncle” phenomenon 
produced by certain schemes such as Roe’s FDS for a supersonic flow over a blunt body, 
first reported by Peery and Imlay [37]. Figure 14 plots the Mach contours and distibution 
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of variables along the symmetry line for a Mach 6 flow. The second-order HUS, using Van 
Leer/Osher splittings, result is free from the “carbuncle” phenomenon and indicates an 
extremely good shock-capturing capability, with nearly no intermediate shock points. The 
convergence history displays monotonic decreasing of residual towards machine zero. 

V.5 Hypersonic Compression-Ramp Problems 

This set of problems consists of hypersonic laminar flows over 15° and 20° wedges. 
The flow conditions for the two test cases are the same. The Reynolds number based on 
the flat plate lenght L is Re L = 18000 and the free stream Mach number = 10. The 
other parameters are: Too = 52 K, T wa u = 290 K , Pr = 0.72, and 7 = 1.4. 

These problems have been extensively studied for they were proposed to the second 
Workshop for Hypersonic Reentry Problems held in Antibes, France in 1991. The Pro- 
ceedings [40] contains several contributions that provide grid-independent solutions for 
comparison. A major difference between these two flows is that the 20° ramp angle pro- 
duces a shock sufficiently strong to cause flow separation of the boundary-layer at the 
ramp. 

These two configurations are calculated on a 175 x 50 grid, see in Fig. 15 the grid 
for 15° ramp. An analysis of the available results [40] leads us to adopt the following 
strategy for mesh generation. The space step in the x-direction is set to be constant with 
Ax/L = 0.01. The vertical spacing is exponentially stretched with a constant parameter 
less than 1.08, the smallest mesh size being A y/L = 0.0002 at the wall. The corresponding 
cell Reynolds number Re c = Re L Ay/L is 3.6. The downstream end of ramp surface 
exhibits an expansion comer with a —20° angle, as proposed by Joulot [25]. This expansion 
comer is intended to simulate as in the experiment the rapid expansion wave taking place 
at the end of the compression ramp. 

Concerning with first the 15° wedge. The Mach number contours are displayed in 
Fig. 15b and the surface quantities in Fig. 16, respectively for (a) pressure coefficient, (b) 
Stanton number, and (c) friction coefficient. Due to the expansion comer, the pressure 
suddenly drops there and the Stanton number first slightly increases and then drops also. 
Included also are results obtained using the Osher-Solomon splitting, nearly in complete 
agreement with the HUS results. 

Next tur nin g to the 20° case. In this case, the increased ramp angle casuses a flow 
separation. In Fig. 17, we compare flow quantities at the wall predicted by the HUS 
method with that by the Osher-Solomon splitting. Both methods again yield virtually 
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identical profiles. This validates that for the boundary-layer type flow, in which the linear 
field dominates, the HUS method recovers as intended the accuracy of the Osher- Solomon 
method. 


VI. Concluding Remarks 

In this paper we presented a new upwind method with detailed mathematical con- 
struction and analysis. The method is based on a field by field decomposition and endows 
the genuinely nonlinear and linearly degenerate waves with effective and compatible pairs 
of FVS and FDS respectively, thus is termed the hybrid upwind splitting (HUS). We be- 
gem by defining the compatibility condition for upwind biasing, resulting in a set of path 
(field)-defined schemes. We showed that the formalism is sufficiently general to include 
existing approximate Riemann solvers, such as the Roe, the Collela-Glaz, and the Osher- 
Solomon methods, and the Van Leer and the Steger- Warming flux vector splittings. We 
gave an interesting interpretation and distinguished the FDS and FVS schemes in terms of 
the path dependency. Further analysis reveals clearly the direct numerical consequences of 
the underlying differences of these schemes, specifically their strengthes and weaknesses. 

Based upon the formalism, we formulated a suitable hybridization that replaced the 
weak link of the basic scheme with a compatible path satisfying such desirable properties 
as: robustness (stability), entropy-satisfying, and low numerical diffusivity. Specifically, 
the genuinely nonlinear subpaths are endowed with a FVS, and the linearly degenerate 
one with a FDS. Explicit formulas using the well-known Osher-Solomon, Collela-Glaz, 
and Roe-type hybridization were given. Unlike the FDS, these expressions do not involve 
evaluations of endpoints of any of the nonlinear subpaths, nor the sonic point. As a result, 
the HUS method is simpler (more efficient) than the underlying FDS, but only adding a 
slight overhead to the FVS in the linearly degenerate subpath, especially as we make use 
of the symmetry property of the FVS with respect to the Mach number. In addition to 
the Van Leer splitting, am one-parameter family of the Boltzmann schemes were derived 
and also included as an alternative of FVS. 

To the Euler and Navier- Stokes equations, we applied the HUS methods on several 
test problems to evaluate their performance against the underlying FDS and FVS methods. 
These tests convincingly revealed the benefits of the proposed HUS methods, in the respects 
of stability for solving a strong shock or rarefaction, accuracy for both inviscid and viscous 
predictions. A formula describing the rate of entropy generation of the numerical schemes 
was given to provide a valuable information about the associated numerical dissipation, 
resulting in an unusual but enlightening insight into a general flux splitting method. 
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The HUS methods in this paper were presented only in the context of explicit schemes. 
However, issues pertaining to linearization must also be considered as an efficient marching 
or iterative scheme to achieve converged solution is preferred. 
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Appendix 1 

For the sake of completness, we apply in this appendix the formal transposition of 
FDS methods in our framework for upwind biasing (proposed in Section II) to some of the 
prominent approximate Riemann solvers : namely the Osher, the Collela-Glaz and finally 
the Roe solvers. The following derivations only intend to recall some of their basic features 
used in Section III. For a complete derivation and furthers results, we refer the reader to 
the original works. 

Al.l The Osher-Solomon Scheme 

This upwind method falls by construction into the formalism proposed in Section II. 
For any given u L and ur in U, $(s;u L ,u R ) is defined to be a Lipshitz continuous path 
made up of p subpaths, the Hh subpath $*(., ui, ur) being a solution of the following 
ODE (see (1.8)) : 

l = rt( $ t{s .,u L ,u R )). (Al.l) 

ds 

Generally speaking, the Cauchy problem (Al.l) admits an unique solution for \u R - u L \ 
sufficiently small. A given ordering of the p subpaths yields the initial condition for each 
of the p ODE above. Let us underline that the resulting path $ heavily depends on the 
ordering under consideration. Two orderings have been singled out in [35], namely the 
natural order $ = u£ =1 $ fc (see also [42]) and the reverse one $ = U \=p$k put forward by 
Osher-Solomon for its theoretical advantages (see Section V for a numerical illustration). 
Actually for the Euler equations and under fairly general thermodynamic assumptions (see 
[35]), both paths exist in the large as long as there is no cavitation. 
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Let us denote |A| : U — * Mat(7i p ,TZ P ) the absolute value of the diagonalizable Jaco- 
bian matrix A. Then, B ± are given by 

B* = |(4 ± |A|), (41.2) 

which obviously satisfy both assumption (2.2) and the compatibility condition (2.3). The 
associated numerical flux therefore reads with clear notation : 

f(u L ,u R ) = ^(f(u L ) + f(«R) - l A *l (A1.3) 

where we have used the following identity deduced from (Al.l) 

(|A| - =0, fors£ [s*,s* +1 [. (A1.4) 

Since all the fields are assumed to be either linearly degenerate or genuinely nonlinear, 
the above vector- valued integration can be explicitly performed and yields 

/(«L,«r) = ^ (f( w L) + f (“ft) - ^2 sigH^kXfOik+i) - 2f(u k+ i) + f(u k ))) , (A1.5a) 

where for A;, 1 < k < p, we have set 

sign(Ak) = lim sign(A k ($(s;uL,UR.))), u k = $(sk;uL,UR.). (A1.56) 

s ^ s k+i 

By convention, u k+ i = u k if A*(ujfc)A*(u*+i) > 0; otherwise u k+ i is the unique point of 
<& k such that A* = 0. 

A1.2 The Collela-Glaz scheme 

In this approach and by contrast with the previous one, the Riemann problem is solved 
in the phase space using only (admissible and nonadmissible) shocks and contact waves. 
The resulting solution is uni- valued (while multi-valued in the Osher- Solomon case, see 
Dubois [11] and Forestier-Le Floch for a precise definition [13]) and is averaged as in the 
Godunov method. For ur and u l close enough, such a solution exists and is made of (at 
most) p waves separated by ( p — 1) constant states ur = u 0 , «i, u p = ur such that 
there exist p real numbers cr k satisfying 

<^jt(^it+i w/fc) = f (ufc+i) f(^fc)? 1 — k — P- (A1.6) 
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For any ul and ur'ihU for which the construction above makes sense, the path $(.; u l,ur) 
is defined to be U^ =1 $jt(.; u L , u R ) where for k, 0 <k<{p- 1) 


with 


*k(s;u L ,u R ) = u t + (7^r7^"‘ +1 ~ Ul) ’ far s e 

_ eL Ww_ z , i<k<p. 

so-0, s t - | BI+ i_ UI i*,’ - 


(A1.7a) 


(A1.76) 


This path is locally Lipshitz continuous. 

We then define for s € [sjt,sjt+i[, 0 < k < p. 


B ± (.,s) = — (cr* ± \(Tk\)IdTit, 


(A1.8) 


which clearly satisfy assumption (2.2). Thanks to the Ranldne-Hugoniot jump relations 
(A1.6), we have 

crk(uk+i — Wfc) = {(ur) — f(«i)i 


Jt=i 

so that (2.3) is valid. The numerical flux can be written as 

p /•**+ 

2 \ ■ -v”—/ z j i 11 o$ 


1/ , . f 9k+l , ,d$k(s-,u L ,u R ) \ 

f(UL,U R ) = 5 (f(u t ) + f(«R) - £ J n & r 

An equivalent form is 

f(u L ,u R ) = x (f( u L) d - f( u ii) sign((Tk)( f ( u k+i) f ( Uk ))) ■ 

^ ^ *-=1 


(A1.9) 


(A1.10) 


Remark 

Al.l- From the two constructions above, the Godunov scheme clearly falls into the 
framework we propose. 

A1.3 The Roe scheme 

For two given states u L and u R in U, the Roe method consists in averaging, as in the 
Godunov method, the exact solution of the following linear problem 


du = 0, t > 0, x e R, 

dt dx 

UL, X < 0 , 
x > 0. 


{j l 


(Al.lla) 
(Al.l 16) 
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Here, f (ul,ur)(u) = B(ul,ur)(u ), for all u in U. The matrix B(ul,ur ) is called a Roe- 
type linearization and is required to satisfy the following conditions for any given ul,ur 
in U 

B(ui,ur) is diagonalizable with real eigenvalues , (A1.12a) 

B{u l , ur)(ur - u L ) = f(u R ) - f(u t ), (A1.126) 

= A(ul) = A(ur). (A1.12c) 

The exact solution shaxes the same structure with that considered by Collela-Glaz. The 
construction we have proposed for that method can therefore be extended in a straight- 
forward way to the Roe method. We just outline this extension. For k, 1 < Jc < p, let 
^k(uL,UR) be the kth. eigenvalue of B(ul,ur) and r*(tt £, ur) the associated right eigen- 
vector. Then by considering the following classical decomposition 

p 

ur- u L = Vkrk{u L , u R ), (A1.13) 

fc=i 

we are led to introduce the following (p — 1) constant states uj, = uo,ui, ..., u p = ur 
by setting for Ar, 1 < k < p, it* = u*_i + ^krk(uL, ur). The path is constructed as in 
(A1.7), (A1.8), B ± are defined to be 

B ± = B(ul,ur ) ± |R|(u L ,Uii)). (A1.14) 

The numerical flux associated with the Roe method can be written 

f(u L ,UR) = ~(f{u L ) + f(tt R ) - j | ^k(uL, «r)I ~~ k ( S *Q s L ' UR '^ ds'j , (A1.15) 

since we have 

(| B\(ul,ur) — \\k(uL,UR)\IdKr ) — q ? for s €[sk,Sk+i[. (A1.16) 

The Roe numerical flux can be rewritten 

f(u L ,UR ) = ]- (f(lt L ) + f(u«) - 53 si S n ( A k(u L , U R ))(f(uL, U R )(u k+1 ) - f(u L , U R )(u k ))] . 

*•=1 / 

(Ai.n) 

Let us stress that the above reconstruction is rather formal. Indeed, the states it* asso- 
ciated with a general Roe method do not a priori necessarily stay in U. In other words, 
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intermediate states are not necessarily kept within their physical range. We stress that 
this difficulty is inherent to the Roe method itself. 

Although the approximate Riemann solution given by the Roe method is only made 
of linear waves, for the sake of convenience, a A'th subpath of Roe s path $ will be said 
to be gen uin ely nonlinear if the associated eigenvalue A *(ul,ur) is consistent with the 
eigenvalue A*(tx) of A{u) such that condition (1.5a) is met. Otherwise if (1.5b) is satisfied, 
then the kth subpath will be referred to as linearly degenerate. 

Rem arks 

A1.2- The analogy between the Osher and Roe flux functions, (A1.5) and (A1.17), 
has been already pointed out by Harten, Lax and Van Leer [22] (see also [36]). 

A1.3- In a recent work devoted to multicomponent Euler equations, Mehlman [33] has 
observed that provided ul and ur are sufficiently close, then the Collela-Glaz method can 
be re-expressed equivalently in terms of a Roe-type linearization method. 

In view of these remarks, we call attention to two main features of the FDS methods 
we have revisited in our framework. First, besides their very similar expressions, their 
formulations are clearly path dependent. Then, when ur and ur are connected in the phase 
space by a single contact discontinuity with zero speed, all the FDS methods considered 
above return the exact solution in the following sense 

«t(“L,“R) = “t, Ur(ul,Ur) = UR (A1.18) 

since we have (see [4], [35], [45]) 

/(ui,ur) = f(ui,) = f(«ft)- (A1.19) 


Appendix 2 

In this Appendix, we derive an one-parameter family of FVS methods, we have referred 
to in Section IV as the T-Boltzmann schemes. Here, T belongs to ]1,3]. These schemes 
come from a set of velocity repartitions parameterized by T and designed by Kamel in [26] 
to be formally consistent with the one-dimensional Euler equation modeling an isentropic 
calorically perfect gas, for which the isentropic coefficient is given by T €]1,3]. These 
velocity functions have bounded support and turn out to fall into the general framework 
due to Perthame [38], [39] for designing entropy- satisfying Boltzmann schemes. Here, we 
apply this framework to the family of velocity distributions quoted above. 


59 



A2.1 Perthame’s kinetic formalism [38] 


This section, is concerned with the notation and a few basic facts about Perthame 
methods. Most of the material here is taken from the work of this author, [38] and [37], 
to which the reader is referred for a comprehensive study and further comments. 

Let 0 : u € H 0(w) € be a nonnegative function satisfying the three consistency 
conditions 

0(— u>) = ©(a;), [ Q(u)du> = 1 , f u 2 Q(uj)du = 1 . (A2.1) 

Jn Jn 

Moreover, we shall assume that the function 0 is boundedly supported, that is, there exists 
ujif (E 7i + such that 

0(u?) = 0, | to |> u ) m- (A2.2) 

In [38], Perthame has shown how to derive from given kinetic entropy pairs such functions 

0 . 


Starting with the given initial data ito(x) = it(0,x), Perthame defines the velocity 
repartition function associated with this data by 

KM) w-KM) 


Ho(x,v) = 


v / KM) 7KM) ^VKMVKM)^ 

, N X , n /KM) ( u> — v(0,x) \_,p(0,x) 

= M0.*)y ^ e(- 7 === / ) - 


with 


A = 


VkMVKM)' 

(3-7) 


(A2.3a) 


A '')"’ (x, u), ( A2.36) 

p(0, x) 


2 ( 7 - 1 )' 

It can be easily checked that the above functions satisfy: 

p(0,x) = / no(x,u>)cL', pw(0, x) — w(0, x) / fjL 0 (x,u>)du, 
Jn Jn 

pv{ 0,x)= / w/io(x,u/)cL;, 

Jn 

pE{ 0,x) = J j^/x 0 (x,o;) + i/ 0 (x,u>)j dw. 

Now, consider the following linear transport problems 

( dtp 4- ud x p = 0, t > 0, x,u> € 72., 

<9*1/ + = 0, 

l /z(0,x,u>) = po(x,u), i/(0, x,u) — uo(x,u:), 


(A2.3c) 

(A2.4a) 

(A2.46) 

(A2.4c) 


(A2.5) 
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whose exact solutions are explicitly known and given by 

n(t,x,u>) =no(x-ut,u), v(t,x,u}) = v 0 (x-ut,u), Vi > 0, *,« € 'll. (A2.6) 

Now let h be a small time step, then 

p{t,x) = f p(t,x,u})d<jj, pw(t,x) = w(t,x) I p,(t,x,oj)duj, ( A2.7a ) 

Jn Jv - 

pv(t,x ) = f u)p(t,x,u>)du>, (A2.7b) 

Jn 

pE{t, x) = /,( Y/i(i,i ) w) + K*»*. w )} <iw * (A2.7c) 

axe first-order in h , approximations of the solution to the Euler equations for t < h (cf. 
[38]). 

Due to the very similarity of the formal construction proposed by Perthame with the 
collisionless Boltzmann equation, the above velocity repartition functions provide us with 
a kinetic interpretation for the Euler equations. 

A Boltzmann scheme is defined to be a scheme in conservation form equipped with 
the following numerical flux 

/(u l ,«r) = / + Ol) + /“(“*). V(u L ,u*)eW (A2.8) 

where the split fluxes (/+/“) axe shown to be given for any given piecewise constant data 
of the form (1.9) by 

f+(u) = j + ( 0 , 0 ,u;)i/o(x,u;)jda;, (A2.9a) 


and r 3 1 

f-(u) = J J(u;,u; 2 ,^-j^o(x,u;) + ( 0 , 0 ,w)i/ o (x,<*>)jdu>, 

provided that the following CFL like-condition is met : 




(A2.96) 


(A2.10) 


The condition stated above clearly enforces each velocity distribution of a given computa- 
tional cell to stay within the adjacent cells. 

Notice that summing (A2.9a) with (A2.9b) gives the consistency condition 

f + (u) + /"(«) = %), Vuew. 
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Boltzmann schemes are shown in [38] to preserve the positivity of the density p and of 
the pressure p under the CFL-like condition (A2.10). Moreover, these schemes are entropy 
satisfying as soon as the generating function 0 is deduced from a given kinetic entropy. 

A2.2 Kaniel velocity distribution functions [26] 


We now turn to describing a particular set of compactly supported velocity distribu- 
tions. Proposed by Kaniel [26], these distributions can be actually defined by introducing 
the following one parameter family of functions (®r) 1<r<3 , ©r : G Tl — ► 0 r (w) € Tl + 

0 K") = (A2.11a) 

where for a > 0, denotes the step function 




1, for all |u;| < a. 
i 

^ 0, otherwise. 


Here, the normalization constant is set to 


(A2.116) 


C r = (A2.11c) 

With this choice of Cr, it is a simple matter to verify that for any given T €]1,3], we have 

[ ©r (cj)du; = 2 f ©r(u = 1, 

Jn Jo 

f 2r\ I \ i 2r^T r±l 

I u> ©p(u;)au; = — — — / u> r - l cLj = 1, 

Jn 1 — 1 Jo 


so that requirements (A2.1) and (A2.2) are satisfied. The one parameter family under 
consideration therefore yields velocity repartition functions in the sense of the Perthame 
kinetic formalism. 


Here in fact, we had to slightly modify the original functions introduced by Kaniel 
to make them fall into the framework described above. Instead of (A2.ll), the generating 
function ©r is defined in [26] for two given positive real numbers p and S as follows 

0rk S](u>) = ^Cr[S]| W |^l [ _ v1f ^^ j+vT ^ ] ( t ,), (,42.12a) 

where the normalization constant Cp[5] is chosen to be 

Cr[S] = (rS)^T. (A2.l2b) 


62 



As a consequence, the resulting (l,w 2 ) moments of 0 r axe now respectively given by p 
and Sp r . 

These results allowed Kaniel to define for a given p and S the following quantity 


p = Sp r , so that Ctt"') = FSp r 1 = — > 0. (A2.13) 

\op/ s p 

It is now clear that with respect to the T-isentropic Euler equations, the generating function 
(A2.ll) can be understood as a thermodynamic function depending solely on the density 
p and the entropy S such that its second moment returns the pressure p = Sp r . For in- 
stance, the choice F — 3 yields a thermodynamic function consistent with the monoatomic 
isentropic one-dimensional gas dynamics, while the choice T = 2 can be formally related 
to the shallow water equations. 


We now derive the explicit flux formulas for the Boltzmann schemes based on the 
generating functions (A2.ll). This is the purpose of the following statement. 


Lemma Let be given T g]1,3]. Let c and M respectively denote an associated sound 
speed and a Mach number given by 



(A2.14) 


Then, under the CFL condition 



(A2.15) 


the generating function 0p ( A.2.11 ) yields a consistent FVS numerical flux defined by the 
following pair of split fluxes (/ + ,/~) • 


if M > +1 / + (u) = f(u), /“(«) = 0 ; 
if M < -1 f + (u ) = 0, /“(«) = f(u) ; 


and when M G [—1,1], 


^W = ±f(f4r ±Af + FTT |M|fS ) 


2 vr + 1 


r + i 

f%{») = v ff( u ) + \ (i ± ffy( 2r ~ ( r _ x ) 

r± f,.\ — ( v i \F\ f ±( 


M 




=(-=-+ i-jfrw 


(A2.16a) 

(A2.16&) 


(A2.16c) 


1 / r 3 m m 2 (5r 2 - sr + 3) m \ 

± 2V3r-i ± 2 + r + r 2r 2(3r-i) 1 1 ’r 
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Proof of the Lemma 

In what follows, we only address the derivation of the forward split flux / + . The 
associated backward split flux is then deduced from the consistency relation (2.35). 

Starting with the first component of the forward split flux / + , definition (A2.9a) yields 



_Cr P_ f 

r - 1 y/p/p Jw>0 


u) — V 3-r 

y/p/p 


1 [-vT,+vT]( 


(jJ — V 

Vp/p 


)du>. 


(A2.17) 


A natural change of variables is given by u = zy/Tp/p + v and leads us to introduce the 
following notation 



(A2.18) 


where c and M respectively stand for the sound speed and the Mach number given in 
(A2.14). The above change of variable clearly shows that if M < —1 then /+ (u) = 0. 
Conversely, if M > 1 then f+(u) = pv. 


Now, focusing our attention on M € [—1, 1], the integral (A2.17) becomes 


4 + (“) t-AS-u 


I d z + M 


f 

J-M 


dz 


}• 


(A2.19) 


The last integral of the right hand-side above is explicitly given by 



r ~ l dz = 


r-i 


(l + sign(M)|M|i*r), 


(A2.20) 


where by convention, sign(M) = +1 if M > 0 and sign(M) = —1 otherwise. 

Concerning the first integral in (A2.19), it can be shown that 

jT * |* M|^). (A2.21) 


Therefore when M € [—1,1], the mass forward flux reads 

f t {u) = f (rfr + M + FTT 1 M |fS ) 

which is the expected result. 

Turning to evaluating the second component of the forward split flux, we have 

* w = -hi/ 1 7m ( " 2 - 22) 
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Applying the previous change of variables yields for M € [—1, 1] 


/?» = TZ 


w: 


z 2 I z |r=T dz + 2M / Z I 2 |rrr 
M •'-** 


$ dz + M 2 [ |*|^ dz,} 

J-M 


where the first integral is computed as 


r-i 

2r 

r-i 


z 2 \z\ ^dz = - (l + sign(M) I M |*=r) 


(l + M|Af |^). 


(A2.23) 


In view of (A2.20) and (A2.21), the momentum numerical flux reduces to the following 
formula when M € [—1,1] '• 

/» = »/;« + K 1 + ffi (2r - 1 < r - x) i M |K) ) p - 

We now conclude the proof in deriving the last component of the forward split flux. 
Definition (A2.9a) provides us with 


foM - r -i 2 1 |r " VvnW 


+ A V -ft(u). 


(A2.24) 


Applying again the same change of variables as above, and after some easy algebraic 
calculations, we can split the first integral of (A2.24) into 


t*m + 2 (£t ){£ 


z 3 | z I 1 - 1, dz+3M / z 2 | z dz+2M 2 

M J~ M 


3-r ^ 

z I z <fzj. 


Here, the integral to be calculated is explicitly given by 


/ 1 * = 11 ^( 1 - 


= ^ l ( 1_Af2|M|fS ) (A2 ' : 

Gathering the relations (A2.21), (A2.23) and (A2.25) yields the required expression^ 


Appendix 3 
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This appendix is devoted to a brief derivation of the formula (5.2) stated in Section 
V, accounting for the rate of entropy dissipation /production of a given 3-point scheme. 
Some co mm ents are given, intending to illustrate its interesting property. 

We focus our attention on the system of conservation laws (1.1) that possesses an 
entropy pair (S', F) : U x U — *• (S(u), F(u)) € 7£ x 71, in the sense of Lax (see [28] for 
instance) : 


S is a convex function of it, i.e. S uu > 0, (A3. la) 

5 and F satisfy the following compatibility condition S u f u = F u . (A3.16) 


In view of the requirement (A3. lb), we see that any given smooth solution u of (1.1) 
satisfies the following additional conservation law 


dS(u) dF (u) 
dt dx 


= 0 . 


(A3. 2) 


Weak solutions of (1.1) no longer satisfy this additional conservation law because the jump 
relation associated with (A3. 2) is not compatible, generally speaking, with the Rankine- 
Hugoniot conditions coming from (1.1), in the sense that it is not satisfied for the same 
value of the discontinuity velocity. This property yields a suitable criterion for selecting 
physically relevant weak solutions. 

Since weak solutions of (1.1) are not uniquely determined by their initial data (1.2), 
the physically relevant solutions axe selected to be the limit, as the viscosity e goes to zero, 
of solutions u ( of the viscous system 


du df (u) dPu 

at + dx ~ e !hP = °’ e > °* 


(A3.3) 


Limit solutions of (A3.3) satisfy in the sense of distributions, the following inequality [20] 


dS(u) + dF («) ^ Q 


(A3 .4a) 


dt dx 

That is, if u is piecewise smooth with discontinuities, then (A3.2) holds pointwise in the 
smooth regions, while across a discontinuity 


-<t{s(uk) - 5(u L )} + {f(ur) - F(u l )} < 0, 


(A3. 46) 


where o denotes the discontinuity velocity. Relations (A3.4b) are called, after Lax, the 
entropy conditions. We note that for a shock wave, namely a discontinuity associated with 
a genuinely nonlinear field, inequality (A3.4b) must be satisfied in a strict manner. 
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When concerned with the numerical approximation of weak solutions of system (1.1), 
we must make sure that the method selects only physically relevant solutions, namely the 
weak solutions that satisfy entropy in the sense of (A3.4). In this way, we ask the method 
(1.9)-(1.14) to be consistent with the entropy condition in the following sense. 

For the method under consideration, there exists a numerical entropy-flux F :UxU — *• 
F(u[,ur ) € F assumed to be Lipschitz-continuous and consistent with the entropy-flux : 

F(u, u) = F(u), Vu € U , (A3. 5a) 

such that the family of approximate solutions u h defined in (1.9) satisfies 

S(u" +1 ) - s«) + 7^+1 - f?_l) < o , iez,neM, (AZ.bb) 

fl x 2 2 

where we have set 

«F(ur+i,u?). (A3.5 c) 

The relevance of requirement (A3.5) is assessed by an extension of the Lax- Wendroff Theo- 
rem [20] which states that limits in the sense of bounded, Lj oc convergence, of approximate 
solutions given by (1.9)-(1.14) satisfying (A3.5b) axe entropy weak solutions of (1.1) and 
(A3.4). 

Rem ark 

A3.1- Inequality (A3. 5b) can be read as the integral entropy law (A3.4a) applied over 
the computational cell. 

Equipped with these definitions, we see that the entropy condition stated in (A3. 5) 
requires the definition of a numerical entropy-flux, whose derivation actually depends on 
the n um erical flux / itself. Then according to the above remark, the numerical entropy 
inequality (A3. 5b) is the sum of the entropy rates induced by the averaging (2.5a) and 
(2.5b) respectively written for the (i + 1/2) and (i — 1/2) interfaces. (A3. 5b) therefore 
accounts for two different jump resolutions. 

Here, we axe merely interested in the entropy rate created in resolving a single jump 
(u t ,u fi ) in order to avoid a possible misleading mixing of entropy production/dissipation 
bv averaging neighboring (but distinct) contributions. Moreover, we want to circumvent 
the nontri vied issue of deriving the entropy numerical flux for a specific flux splitting method 
and instead we intend to arrive at a unique formula as free as possible of the underlying 
method. These two goals can be achieved in the following way. The derivation below is 
formal but can be nevertheless rigorously justified on the basis of (A3. 5). 
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We first foe vis our attention on the half computational cell (ih x ,(i+l/2)h x )x ( nh , (n+ 
l)h). Applying formally the integral conservation law (1.1) over this domain yields under' 
the CFL condition (2.4), the averaging (5.1a) : 

= <- £(/(«•*>«?+: l) - fW))- (^-«) 

Ag ain formally, in view of Remark 1, applying now the integral entropy law (A3.4a) over 
the same domain gives under the CFL condition (2.4) 


Using the consistency condition (A3.5a), we clearly get : 

$(«*.(<*?. »«■•)) - $(«?) + j£(F(u?,< +1 ) - F(u?)) < 0. 

A similar derivation can be performed also for the rectangle ((i + 1/2 )h x ,(i + 
(nh,(n + l)h), providing us successively with 


(A3. 7a) 


(A3. 76) 
l)h x ) x 


and 


2h 




2h 


S(Sr«, < + 1 ) - S« +1 ) + r (F(u? +I ) - F(“".“r+i)) < 0. 


(A3. 8) 
(A3.9) 


by using the same arguments as above. 

It suffices to stun (A3. 7) and (A3. 9) to deduce the expected relation (5.2) for the 
entropy source term 


r(S(«R(u”,u” +1 )) + S(tt£(u”,u” + 1)) — ^(wf+i) — S( u ?)) 


H (F(u” +1 ) — F(u”)) = Entropy source term. (A3.10) 

h x 


Relation (A3. 10) clearly accounts for the rate of entropy created in the shifted rectangle 
(( ih x ,(i + 1)6*)) x (0, h) in resolving the single jump (it”, u" +1 ). Since such a rate is closely 
related to the amount of numerical dissipation involved in the resolution of that jump (see 
[5], [6]), formula (A3. 10) provides us with a valuable tool, delivering a clear insight into a 
general flux splitting method. 
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Figure 2. Shock tube problem, Case 1: (a) Comparison of the results from the methods of 
Roe, Osher (natural ordering), HUS (Van Leer-Osher hybridization), HUS (Boltzmann- 
Osher hybridization with T = 3), and the exact Riemann solver (solid line). 
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Figure 2. (continued) (b) Comparison of results from the methods of HUS (Boltzmann- 
Osher hybridization with T = 2), Boltzmann (r = 2), Van Leer and the exact Riemann 
solver (solid line). 
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Figure 3. Shock tube problem, Case 2: (a) Comparison of results from the methods of 
Roe, Osher (natural ordering), HUS (Van Leer-Osher hybridization), HUS (Boltzmann- 
Osher hybridization with T = 3) and the exact Riemann solver (solid line). The entropy 
source term is defined in (5.2) and the entropy dissipation is the time history of this term 
measured at x = 0.5. 
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Figure 3. (continued) (b) Comparison of results from the methods of HUS (Boltzmann- 
Osher hybridization with P = 3), Boltzmann (r = 2), Van Leer and the exact Riemann 
solver (solid line). 
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Figure 4. Shock tube problem, Case 3: (a) Comparison of results from the methods of Roe, 
Osher (natural ordering), HUS (Van Leer-Osher hybridization), HUS (Boltzmann-Osher 
hybridization with T = 3), and the exact Riemann solver (solid line). 
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Figure 4. (continued) (b) Comparison of results from the methods of HUS (Boltzmann- 
Osher hybridization with T = 2), Boltzmann (r = 3), Boltzmann (r = 2), Van Leer and 
the exact Riemann solver (solid line). 
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Figure 5. Shock tube problem, Case 4: (a) Comparison of results from the methods of 
Roe, Osher (natural ordering), HUS (Van Leer-Osher hybridization), Van Leer, and the 
exact Riemann solver (solid line). 
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Figure 6. Shock tube problem, Case 5: (a) Comparison of results from the methods of 
Roe, HUS (Van Leer-Osher hybridization with natural ordering), HUS (Van Leer-Osher 
hybridization with reverse ordering), Van Leer, and the exact Riemann solver (solid line). 
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Figure 7. Quasi one-dimensional hypersonic conical flow problem, Moo - 7.95, 
Re = 4.2 x 10 5 , and half cone angle=10°. 
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Figure 8. Supersonic shock wave/laminar boundary layer interaction problem: surface 
pressure, skin-friction coefficient, and convergence history. 
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friction, and. heat transfer by the Van Leer method. 




Figure 11. Hypersonic hyperboloid problem: grid refinement study of surface pressure, 
friction, and heat transfer by the HUS (VL-O) method. 
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Figure 12. Hypersonic hyperboloid problem: comparison of the Van Leer, Osher, and 
HUS solutions on surface pressure, friction, and heat transfer on two fine grids. 
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Figure 14. Inviscid HUS solution of M 6 flow over a circular body. 


86 



Figure 15. Hypersonic laminar flow over a 15' ramp : (a) 174 x 49 Grid (skipped every 
other lines for clarity), and (b) Mach number contours . 
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